Variational Schemes and Geometric Simulations for a Hydrodynamic-Electrodynamic Model of Surface Plasmon Polaritons
Qiang Chen, Lifei Geng, Xiang Chen, Xiaojun Hao, Chuanchuan Wang, Xiaoyang Wang
AAPS
Variational Schemes and Geometric Simulations for aHydrodynamic-Electrodynamic Model of Surface PlasmonPolaritons
Qiang Chen ∗ State Key Laboratory of Complex ElectromagneticEnvironment Effects on Electronics and Information System,Luoyang, Henan 471000, China andUniversity of Science and Technology of China, Hefei, Anhui 230026, China
Lifei Geng, Xiang Chen, Xiaojun Hao, Chuanchuan Wang, and Xiaoyang Wang
State Key Laboratory of Complex Electromagnetic Environment Effectson Electronics and Information System, Luoyang, Henan 471000, China (Dated: February 21, 2019)
Abstract
A class of variational schemes for the hydrodynamic-electrodynamic model of lossless free elec-tron gas in a quasi-neutral background is developed for high-quality simulations of surface plasmonpolaritons. The Lagrangian density of lossless free electron gas with a self-consistent electromag-netic field is established, and the dynamical equations with the associated constraints are obtainedvia a variational principle. Based on discrete exterior calculus, the action functional of this systemis discretized and minimized to obtain the discrete dynamics. Newton-Raphson iteration and thebiconjugate gradient stabilized method are equipped as a hybrid nonlinear-linear algebraic solver.Instead of discretizing the partial differential equations, the variational schemes have better numer-ical properties in secular simulations, as they preserve the discrete Lagrangian symplectic structure,gauge symmetry, and general energy-momentum density. Two numerical experiments were per-formed. The numerical results reproduce characteristic dispersion relations of bulk plasmons andsurface plasmon polaritons, and the numerical errors of conserved quantities in all experiments arebounded by a small value after long term simulations.
PACS number(s): ∗ [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] F e b . INTRODUCTION In the past two decades, there have been impressive developments and significant ad-vancement in applications of Surface Plasmon Polaritons (SPPs), bringing many new ideasinto traditional electromagnetics and optics, such as the lithography beyond the diffractionlimit, chip-scale photonic circuits, plasmonic metasurfaces, bio-photonics, etc. [1–10]. Inthe field of metal optics, plasmonics focuses on the collective motions of free electron gasin metal with self-consistent and external electromagnetic fields whose first-principle modelis the classical particle-field theory [1–3]. Direct applications of the first-principle modelin 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 plasmonicphenomena [2, 3]. In a mesoscopic context, kinetic and hydrodynamic descriptions are basicphysical models of plasmonics, which involve both the dynamics of free electron gas and anelectromagnetic field. Therefore, high-quality numerical schemes and simulations based onthe hydrodynamic model are necessary in plasmonic research.The physics of SPPs can be described by the free electron gas model whose dynamicalequations are hydrodynamic and Maxwell’s equations [1]. For Maxwell’s equations, manynumerical schemes, such as the Finite-Difference Time-Domain (FDTD) method, the FiniteElement (FE) method, and the Method of Moments (MoM) have been developed, that arewidely used in modern electromagnetic engineering, Radio Frequency (RF) and microwaveengineering, terahertz engineering, optical engineering, metamaterial design, accelerator de-sign, fusion engineering, radio astrophysics, geophysics, and even biomedicine [11–14]. Forthe hydrodynamic equations, there are also many well-designed numerical schemes. In ad-dition to the Finite Difference (FD) method and the FE method, the Finite Volume (FV)method for conservation-type equations is very popular in Computational Fluid Dynamics(CFD), which is expected to achieve superior conservations, e.g. total energy conservation[15]. Because of their nature and the advantages in their respective fields, traditional nu-merical methods are often extended and reformed to simulate more complicated or hybridsystems [10, 16]. When it comes to plasmonic phenomena, the simplest numerical treatmentis introducing the DL model, which is a local linear response of the hydrodynamic equationsinto the FDTD or FE methods [2, 3]. This numerical model can be conveniently solved by cir-2ular convolution or Auxiliary Differential Equation (ADE) techniques [10]. This linearizeddispersion model provides us with abundant information about plasmonic perturbations,such as the linear dispersion relations and polarization modes, but the important nonlinearand dynamical properties, such as mode mixing, High Harmonic Generation (HHG), andthe Kerr effect, are ignored [6, 9]. An FDTD-type method for Hydrodynamic-Maxwell equa-tions of plasmonic metasurfaces can also be found, and many interesting results about thenonlinear effect are shown in the simulations [10]. Previous work showed that plasmonicphenomena are attractive and important, which means that further numerical research isnecessary. As a basic tool for computational plasmonics, more advanced numerical schemesfor the free electron gas model are needed.Because of the nonlinearity and the multi-scale nature of the Hydrodynamic-Maxwellequations, high-quality simulations of SPPs face challenges. For example, numerical er-rors involving the momentum and energy of the electron gas and the electromagnetic fieldcan coherently accumulate, though these errors may be small in each numerical step. Thebreakdown of conservation laws over a long simulation time amounts to pseudo physics. It isdesirable, therefore, to use numerical integrators with good global conservative properties.The canonical symplectic integrators for Hamiltonian systems with a canonical structurefirst developed by K. Feng et. al. are known as a class of structure-preserving geometricalgorithms that have excellent numerical performance in long-term simulations [17–28]. Un-fortunately, the hydrodynamic system does not possess a simple canonical structure, whichmeans the canonical symplectic integrators do not apply directly [29]. As an alternativemethod, J. Marsden and M. West developed the variational integrator based on the discreteHamiltonian principle for Lagrangian systems [30–32]. The discrete variational principlepreserves the Lagrangian structures of dynamical systems in Lagrangian form. As a signifi-cant advantage, the numerical errors of conserved quantities are bounded by a small valueover a long time. The variational integrator has been widely used in many complex systems,especially in geophysics and plasma physics [33–43].In this work, we construct a class of variational schemes for geometric simulations ofSPPs, and we show the advantages of the algorithms via several numerical experiments. InSec.II, the Lagrangian density of lossless free electron gas with a self-consistent electromag-netic field is established, and the dynamics with constraints are obtained via the variationalprinciple. With an appropriate derivation, the standard Hydrodynamic-Maxwell equations3re obtained in an arbitrary gauge. In Sec.III, the action functional is discretized via Dis-crete Exterior Calculus (DEC) and minimized to obtain discrete dynamics. Equipped withefficient nonlinear and linear algebraic solvers, i.e., Newton-Raphson (N-R) iteration andthe Biconjugate-Gradient-Stabilized (BICGSTAB) method, the discrete system is solvable.In App.A, the detailed Jacobian of nonlinear equations is derived, which guides the pro-cedures of nonlinear iteration and associated linear updating. In Sec.IV, the schemes areimplemented to simulate bulk plasmons and SPPs. The numerical results reproduce thecharacteristic dispersion relations of the plasmonic systems, and they exhibit good long-term stability. These desirable features make the algorithm a powerful tool in the study ofSPPs using the hydrodynamic-electrodynamic model.
II. MODEL AND THEORYA. Plasmonics
SPPs are a type of infrared or visible electromagnetic surface waves that travel alonga metal-dielectric interface. The polariton means the wave consists of electron collectivemotion and a self-consistent electromagnetic field between the metal and the dielectric [2, 3].The simplest SPP configuration is shown in Fig.1, where a Transverse Magnetic (TM) modepropagates along an infinite interface between the metal (lower half-space) and the dielectric(upper half-space).
FIG. 1. SPPs configuration at an infinite metal-dielectric interface.
We assume that the media are nonmagnetic. At the region where the electromagneticwave frequency ω is lower than the metal plasma frequency ω p , the metallic permittivity4 < (cid:15) >
0. Then we can write the linearized SPPs fieldby using the Amp`ere-Maxwell equation as, H y ( z ) = ˜ H e ik x x e k z z , (1) E x ( z ) = − i ˜ H k z ω(cid:15) (cid:15) e ik x x e k z z , (2) E z ( z ) = − ˜ H k x ω(cid:15) (cid:15) e ik x x e k z z , (3)for the metal region, and H y ( z ) = ˜ H e ik x x e − k z z , (4) E x ( z ) = i ˜ H k z ω(cid:15) (cid:15) e ik x x e − k z z , (5) E z ( z ) = − ˜ H k x ω(cid:15) (cid:15) e ik x x e − k z z , (6)for the dielectric region, where (cid:15) is the vacuum permittivity, and ˜ H and ˜ H are the magneticfield amplitudes. Applying the interface conditions E x (0) = E x (0) and H y (0) = H y (0)into Eqs. (1)-(6), we can obtain ˜ H = ˜ H and (cid:15) k z = − (cid:15) k z . Then the linearized SPPsdispersion relations are given by using Faraday’s equation as [2], k z = k x − k (cid:15) , (7) k z = k x − k (cid:15) , (8) k x = k (cid:114) (cid:15) (cid:15) (cid:15) + (cid:15) , (9)where k = ω/c is the vacuum wave number.The previous descriptions for SPPs provide us with a basic physical image and importantspectral information. At the low-frequency (terahertz or millimeter) region, Eq. (9) reducesto k x ≈ k and the dielectric confinement strength k z has a very small value, which meansthat the SPPs reduce to weak-bound Sommerfeld-Zenneck (SZ) waves. At the high-frequency(visible and near infrared) region, Eq. (9) reduces to ω ≈ ω sp = ω p / √ (cid:15) , where ω sp is theSurface Plasmon Resonance (SPR) frequency. This means that the SPPs at the resonanceregion are far from the light cone, which has infinitesimal group velocity v g and phasevelocity v p . Both the dielectric and metal confinement strengths k z , of these quasi-staticmodes near the resonance region have very large values, which means the SPPs energy isstrongly localized. As the cornerstone of many plasmonic techniques, the strong-boundsub-wave-length field is an important feature of the SPPs.5 . Lagrangian Form of Lossless Free Electron Gas The Lagrangian density of lossless free electron gas in a quasi-neutral background with aself-consistent electromagnetic field can be defined from the hydrodynamic form as [44, 45], L = L EG + L EM + L Int , (10) L EG = 12 mn v + α (cid:20) ∂∂t n + (cid:53) · ( n v ) (cid:21) − λ (cid:18) ∂∂t µ + v · (cid:53) µ (cid:19) , (11) L EM = (cid:15) (cid:18) − (cid:53) φ − ∂∂t A (cid:19) − µ ( (cid:53) × A ) , (12) L Int = en v · A − e ( n − n ) φ, (13)where the subscripts denote Electron Gas (EG), ElectroMagnetic (EM), and Interaction(Int), respectively. n is the electron density, n is the background particle density, v is theelectron gas velocity, e is the electron charge and m is the electron mass. The electromagneticfield components ( φ, A ) are defined in an arbitrary gauge. This Lagrangian density is givenin constraint form, where the scalar fields α and λ are Lagrangian multipliers. In the lastterm of Eq. (10), Lin’s constraint factor µ is used to establish a complete description for thevelocity and helicity [44, 45].With the Lagrangian density (10), we can construct the action functional S = (cid:82) T (cid:82) Ω L d x d t which involves the physics of lossless free electron gas. With tedious variational calculation,the total variation of S with fixed boundary can be given as, δS = (cid:90) T (cid:90) Ω (cid:26)(cid:18) m v + e v · A − eφ − ∂∂t α − v · (cid:53) α (cid:19) δn + ( mn v + en A − n (cid:53) α − λ (cid:53) µ ) · δ v + (cid:20) en v − (cid:15) (cid:18) (cid:53) ∂∂t φ + ∂ ∂t A (cid:19) − µ (cid:53) × (cid:53) × A (cid:21) · δ A − (cid:20) e ( n − n ) − (cid:15) (cid:18) (cid:53) φ + (cid:53) · ∂∂t A (cid:19)(cid:21) δφ + (cid:20) ∂∂t n + (cid:53) · ( n v ) (cid:21) δα − (cid:18) ∂∂t µ + v · (cid:53) µ (cid:19) δλ + (cid:20) ∂∂t λ + (cid:53) · ( λ v ) (cid:21) δµ (cid:27) d x d t. (14)6aking variational derivative of the action functional S with respect to the fields q =( n, v , A , φ, α, λ, µ ) and minimizing S by leading variational derivatives equal to 0, we canobtain the Euler-Lagrange equations of lossless free electron gas, δSδn = 0 ⇒ ∂∂t α = 12 m v − v · (cid:53) α + e v · A − eφ, (15) δSδ v = 0 ⇒ mn v + en A = n (cid:53) α + λ (cid:53) µ, (16) δSδ A = 0 ⇒ ∂ ∂t A = − (cid:15) µ (cid:53) × (cid:53) × A − (cid:53) ∂∂t φ + e(cid:15) n v , (17) δSδφ = 0 ⇒(cid:53) φ + (cid:53) · ∂∂t A = − e(cid:15) ( n − n ) , (18) δSδα = 0 ⇒ ∂∂t n = − (cid:53) · ( n v ) , (19) δSδλ = 0 ⇒ ∂∂t µ = − v · (cid:53) µ, (20) δSδµ = 0 ⇒ ∂∂t λ = − (cid:53) · ( λ v ) . (21)Eqs. (15) and (21) are Lagrangian multiplier equations. Eq. (16) defines the canonical mo-mentum density of free electron gas. Eqs. (17) and (18) are Maxwell’s equations. Eq. (19) isthe continuity equation of free electron gas. Eq. (20) defines the dynamics of Lin’s constraintfield. These equations provide us with a complete model to describe plasmonic phenomena. C. Hydrodynamic-Maxwell Model
The Euler-Lagrange equations (15)-(21) can be recognized as a general form of the well-known hydrodynamic-Maxwell equations. Here we give a detailed derivation to illustratethe relation between two kinds of equations.7dopting the temporal gauge φ = 0 explicitly, we can define the electromagnetic fieldcomponents as E = − ∂ A /∂t and µ H = (cid:53) × A . Taking the derivative of Eq. (16) withrespect to t , we obtain, m (cid:18) ∂∂t n (cid:19) v + mn ∂∂t v + e (cid:18) ∂∂t n (cid:19) A + en ∂∂t A = (cid:18) ∂∂t n (cid:19) (cid:53) α + n (cid:53) ∂∂t α + (cid:18) ∂∂t λ (cid:19) (cid:53) µ + λ (cid:53) ∂∂t µ. (22)With tedious algebraic calculation, we obtain the momentum equation of free electron gasby substituting Eqs. (15) and (19)-(21) into Eq. (22) as, mn ∂∂t v = ( mn v + en A − n (cid:53) α − λ (cid:53) µ ) · (cid:53) v + ( mn v + en A − n (cid:53) α − λ (cid:53) µ ) × (cid:53) × v + v · (cid:18) λn (cid:53) n (cid:53) µ − (cid:53) λ (cid:53) µ − λ (cid:53) (cid:53) µ − n (cid:53) (cid:53) α + en (cid:53) A (cid:19) + en (cid:20) − ∂∂t A + v × ( (cid:53) × A ) (cid:21) . (23)Then taking the gradient (cid:53) on Eq. (16), we obtain, − mn (cid:53) v = en (cid:53) A − (cid:53) λ (cid:53) µ + λn (cid:53) n (cid:53) µ − λ (cid:53) (cid:53) µ − n (cid:53) (cid:53) α. (24)At last, we obtain the velocity equation of free electron gas by substituting Eqs. (16) and(24) into Eq. (23) as, ∂∂t v = − v · (cid:53) v + em (cid:20) − ∂∂t A + v × ( (cid:53) × A ) (cid:21) . (25)Based on the previous derivation, Eqs. (17)-(19) and (25) make up the standard hydrodynamic-Maxwell equations, which can be rewritten in a compact form as, ∂∂t n + (cid:53) · ( n v ) = 0 , (26) ∂∂t v + v · (cid:53) v = em ( E + µ v × H ) , (27)1 c ∂ ∂t A + (cid:53) × (cid:53) × A = µ en v , (28) (cid:53) · ∂∂t A = − e(cid:15) ( n − n ) . (29)The hydrodynamic-Maxwell equations (26)-(29) are complete and self-consistent. Theelectromagnetic field components enter the hydrodynamic velocity equation (27) via Lorentzforce. The hydrodynamic fields enter Maxwell’s equations (28)-(29) via the current and8harge densities. In these equations, both the dynamics of electron collective motion andself-consistent field are involved. Based on previous discussion, both the hydrodynamic-Maxwell equations (26)-(29) and Euler-Lagrange equations (15)-(21) of lossless free electrongas can be used as the basic physical model of plasmonics. III. NUMERICAL STRATEGIESA. DEC Based Discretization
The traditional numerical methods for an infinite dimensional dynamical system focus onthe solving techniques for relevant differential or integral equations, such as the FDTD, FEMor FV types of numerical schemes for the hydrodynamic-Maxwell equations (26)-(29) or eventhe alternative Eqs. (15)-(21) in plasmonic research. Here we construct a numerical strategyfor hydrodynamic-electrodynamic model based plasmonic phenomena simulations in a differ-ent way. Instead of directly discretizing the previous differential equations, we reconstructthe discrete dynamics via discrete variational principle, which is an infinite dimensionalHamilton’s principle analog on discrete space-time manifold. These variational schemeshave particular advantages in high-quality long term simulations, as they can preserve thediscrete Lagrangian symplectic structure, gauge symmetry and general energy-momentumdensity, which will be discussed in detail in the following sections.The first step in numerical simulation is discretization. As a differential geometry basednumerical framework, DEC defines complete operational rules on a discrete differential man-ifold, which form a cochain complex [42, 46–50]. To solve the continuous system numerically,the space-time manifold is discretized using a rectangular lattice (other lattices are also vi-able). Then the scalar fields, which are 3-forms, i.e. ndx ∧ dy ∧ dz and λdx ∧ dy ∧ dz , naturallylive on the volume center of the discrete spacelike submanifold, n ni + ,j + ,k + : n (cid:18) t n , x i + ∆ x , y j + ∆ y , z k + ∆ z (cid:19) , (30)where ( t n , x i , y j , z k ) is the coordinate of the lattice vertex. The velocity, multiplier, constraintand gauge 1-forms v = v i dx i , αdt , µdt and A = A ν dx ν naturally live along the edges of the9pace-time lattice, φ n + i,j,k : φ (cid:18) t n + ∆ t , x i , y j , z k (cid:19) , (31) A nxi + ,j,k : A x (cid:18) t n , x i + ∆ x , y j , z k (cid:19) , (32) A nyi,j + ,k : A y (cid:18) t n , x i , y j + ∆ y , z k (cid:19) , (33) A nzi,j,k + : A z (cid:18) t n , x i , y j , z k + ∆ z (cid:19) . (34)In the above discretization, a half integer index indicates along which edge does the fieldresides. Then the 2-forms, e.g., F = d A , and the 4-forms, e.g., F ∧ ∗ F , are defined onthe faces and volume center of the space-time lattice respectively, where d is the exteriorderivative operator. The discretization of the spacelike submanifold is shown in Fig.2. FIG. 2. DEC-based discretization of the hydrodynamic-electrodynamic model in a rectangularlattice. The discrete spacelike submanifold is shown, and other types of discrete submanifolds canbe given in the same way.
Based on these definitions, the exterior derivatives of discrete differential forms are natu-rally obtained via the forward difference operators. The exterior derivatives on the spacelikesubmanifold are given as, d φ = ( (cid:53) φ ) i dx i = φ i +1 ,j,k − φ i,j,k ∆ x dx + φ i,j +1 ,k − φ i,j,k ∆ y dy + φ i,j,k +1 − φ i,j,k ∆ z dz, (35)10 A = ( (cid:53) × A ) i dx j ∧ dx k = (cid:32) A zi,j +1 ,k + − A zi,j,k + ∆ y − A yi,j + ,k +1 − A yi,j + ,k ∆ z (cid:33) dy ∧ dz + (cid:32) A xi + ,j,k +1 − A xi + ,j,k ∆ z − A zi +1 ,j,k + − A zi,j,k + ∆ x (cid:33) dz ∧ dx + (cid:32) A yi +1 ,j + ,k − A yi,j + ,k ∆ x − A xi + ,j +1 ,k − A xi + ,j,k ∆ y (cid:33) dx ∧ dy, (36) d ∗ A = (cid:53) · A dx i ∧ dx j ∧ dx k = (cid:32) A xi + ,j,k − A xi − ,j,k ∆ x + A yi,j + ,k − A yi,j − ,k ∆ y + A zi,j,k + − A zi,j,k − ∆ z (cid:33) dx ∧ dy ∧ dz, (37)where ∗ is the Hodge star operator, which generates the Hodge dual form of the primarydiscrete form.By using the DEC, the Lagrangian density (10) in the interval [ t n , t n +1 ] can be discretizedas, L n ∼ n +1 di,j,k = L n ∼ n +1 d EG i,j,k + L n ∼ n +1 d EM i,j,k + L n ∼ n +1 d Int i,j,k , (38) L n ∼ n +1 d EG i,j,k = 12 mn ni + ,j + ,k + (cid:16) v n xi + ,j,k + v n yi,j + ,k + v n zi,j,k + (cid:17) + α n + i,j,k (cid:32) n n +1 i + ,j + ,k + − n ni + ,j + ,k + ∆ t + n n +1 i + ,j + ,k + v n +1 xi + ,j,k − n n +1 i − ,j + ,k + v n +1 xi − ,j,k ∆ x + n n +1 i + ,j + ,k + v n +1 yi,j + ,k − n n +1 i + ,j − ,k + v n +1 yi,j − ,k ∆ y + n n +1 i + ,j + ,k + v n +1 zi,j,k + − n n +1 i + ,j + ,k − v n +1 zi,j,k − ∆ z (cid:33) − λ ni + ,j + ,k + µ n + i,j,k − µ n − i,j,k ∆ t + v nxi + ,j,k µ n − i +1 ,j,k − µ n − i,j,k ∆ x + v nyi,j + ,k µ n − i,j +1 ,k − µ n − i,j,k ∆ y + v nzi,j,k + µ n − i,j,k +1 − µ n − i,j,k ∆ z , (39)11 n ∼ n +1 d EM i,j,k = (cid:15) φ n + i +1 ,j,k − φ n + i,j,k ∆ x + A n +1 xi + ,j,k − A nxi + ,j,k ∆ t + φ n + i,j +1 ,k − φ n + i,j,k ∆ y + A n +1 yi,j + ,k − A nyi,j + ,k ∆ t + φ n + i,j,k +1 − φ n + i,j,k ∆ z + A n +1 zi,j,k + − A nzi,j,k + ∆ t − µ (cid:32) A nzi,j +1 ,k + − A nzi,j,k + ∆ y − A nyi,j + ,k +1 − A nyi,j + ,k ∆ z (cid:33) + (cid:32) A nxi + ,j,k +1 − A nxi + ,j,k ∆ z − A nzi +1 ,j,k + − A nzi,j,k + ∆ x (cid:33) + (cid:32) A nyi +1 ,j + ,k − A nyi,j + ,k ∆ x − A nxi + ,j +1 ,k − A nxi + ,j,k ∆ y (cid:33) , (40) L n ∼ n +1 d Int i,j,k = en ni + ,j + ,k + (cid:16) v nxi + ,j,k A nxi + ,j,k + v nyi,j + ,k A nyi,j + ,k + v nzi,j,k + A nzi,j,k + (cid:17) − e (cid:16) n ni + ,j + ,k + − n (cid:17) φ n + i,j,k . (41)Then the action functional of discrete dynamical system is, S d = N (cid:88) n =0 L n ∼ n +1 d ∆ t, (42) L n ∼ n +1 d = (cid:88) i,j,k L n ∼ n +1 di,j,k ∆ x ∆ y ∆ z, (43)where the functional L n ∼ n +1 d is the Lagrangian of discrete dynamical system in interval[ t n , t n +1 ]. B. Variational Schemes
With DEC-based discretization, we reconstruct a field theory on the discrete space-timemanifold. Then the variational derivatives of action S d with respect to fields reduce topartial derivatives with respect to discrete differential forms. By minimizing the action, we12btain the discrete dynamical equations, ∂S d ∂n ni + ,j + ,k + = ∂ ( L n − ∼ nd + L n ∼ n +1 d ) ∂n ni + ,j + ,k + ∆ t = 0 ⇒ α n + i,j,k − m ∆ t (cid:16) v n xi + ,j,k + v n yi,j + ,k + v n zi,j,k + (cid:17) − e ∆ t (cid:16) v nxi + ,j,k A nxi + ,j,k + v nyi,j + ,k A nyi,j + ,k + v nzi,j,k + A nzi,j,k + (cid:17) + e ∆ tφ n + i,j,k + ∆ t ∆ x v nxi + ,j,k (cid:16) α n − i +1 ,j,k − α n − i,j,k (cid:17) + ∆ t ∆ y v nyi,j + ,k (cid:16) α n − i,j +1 ,k − α n − i,j,k (cid:17) + ∆ t ∆ z v nzi,j,k + (cid:16) α n − i,j,k +1 − α n − i,j,k (cid:17) = α n − i,j,k , (44) ∂S d ∂v nxi + ,j,k = ∂ ( L n − ∼ nd + L n ∼ n +1 d ) ∂v nxi + ,j,k ∆ t = 0 ⇒ mn ni + ,j + ,k + v nxi + ,j,k + en ni + ,j + ,k + A nxi + ,j,k = n ni + ,j + ,k + α n − i +1 ,j,k − α n − i,j,k ∆ x + λ ni + ,j + ,k + µ n − i +1 ,j,k − µ n − i,j,k ∆ x , (45) ∂S d ∂v nyi,j + ,k = ∂ ( L n − ∼ nd + L n ∼ n +1 d ) ∂v nyi,j + ,k ∆ t = 0 ⇒ mn ni + ,j + ,k + v nyi,j + ,k + en ni + ,j + ,k + A nyi,j + ,k = n ni + ,j + ,k + α n − i,j +1 ,k − α n − i,j,k ∆ y + λ ni + ,j + ,k + µ n − i,j +1 ,k − µ n − i,j,k ∆ y , (46) ∂S d ∂v nzi,j,k + = ∂ ( L n − ∼ nd + L n ∼ n +1 d ) ∂v nzi,j,k + ∆ t = 0 ⇒ mn ni + ,j + ,k + v nzi,j,k + + en ni + ,j + ,k + A nzi,j,k + = n ni + ,j + ,k + α n − i,j,k +1 − α n − i,j,k ∆ z + λ ni + ,j + ,k + µ n − i,j,k +1 − µ n − i,j,k ∆ z , (47)13 S d ∂A nxi + ,j,k = ∂ ( L n − ∼ nd + L n ∼ n +1 d ) ∂A nxi + ,j,k ∆ t = 0 ⇒ (cid:15) ∆ t A n +1 xi + ,j,k = − (cid:15) ∆ t ∆ x (cid:16) φ n + i +1 ,j,k − φ n + i,j,k − φ n − i +1 ,j,k + φ n − i,j,k (cid:17) + A nxi + ,j,k +1 − A nxi + ,j,k + A nxi + ,j,k − µ ∆ z + A nxi + ,j +1 ,k − A nxi + ,j,k + A nxi + ,j − ,k µ ∆ y − A nzi +1 ,j,k + − A nzi +1 ,j,k − − A nzi,j,k + + A nzi,j,k − µ ∆ x ∆ z − A nyi +1 ,j + ,k − A nyi +1 ,j − ,k − A nyi,j + ,k + A nyi,j − ,k µ ∆ x ∆ y + en ni + ,j + ,k + v nxi + ,j,k + 2 (cid:15) ∆ t A nxi + ,j,k − (cid:15) ∆ t A n − xi + ,j,k , (48) ∂S d ∂A nyi,j + ,k = ∂ ( L n − ∼ nd + L n ∼ n +1 d ) ∂A nyi,j + ,k ∆ t = 0 ⇒ (cid:15) ∆ t A n +1 yi,j + ,k = − (cid:15) ∆ t ∆ y (cid:16) φ n + i,j +1 ,k − φ n + i,j,k − φ n − i,j +1 ,k + φ n − i,j,k (cid:17) + A nyi +1 ,j + ,k − A nyi,j + ,k + A nyi − ,j + ,k µ ∆ x + A nyi,j + ,k +1 − A nyi,j + ,k + A nyi,j + ,k − µ ∆ z − A nxi + ,j +1 ,k − A nxi − ,j +1 ,k − A nxi + ,j,k + A nxi − ,j,k µ ∆ y ∆ x − A nzi,j +1 ,k + − A nzi,j +1 ,k − − A nzi,j,k + + A nzi,j,k − µ ∆ y ∆ z + en ni + ,j + ,k + v nyi,j + ,k + 2 (cid:15) ∆ t A nyi,j + ,k − (cid:15) ∆ t A n − yi,j + ,k , (49) ∂S d ∂A nzi,j,k + = ∂ ( L n − ∼ nd + L n ∼ n +1 d ) ∂A nzi,j,k + = 0 ⇒ (cid:15) ∆ t A n +1 zi,j,k + = − (cid:15) ∆ t ∆ z (cid:16) φ n + i,j,k +1 − φ n + i,j,k − φ n − i,j,k +1 + φ n − i,j,k (cid:17) + A nzi,j +1 ,k + − A nzi,j,k + + A nzi,j − ,k + µ ∆ y + A nzi +1 ,j,k + − A nzi,j,k + + A nzi − ,j,k + µ ∆ x − A nyi,j + ,k +1 − A nyi,j − ,k +1 − A nyi,j + ,k + A nyi,j − ,k µ ∆ z ∆ y − A nxi + ,j,k +1 − A nxi − ,j,k +1 − A nxi + ,j,k + A nxi − ,j,k µ ∆ z ∆ x + en ni + ,j + ,k + v nzi,j,k + + 2 (cid:15) ∆ t A nzi,j,k + − (cid:15) ∆ t A n − zi,j,k + , (50)14 S d ∂φ n + i,j,k = ∂L n ∼ n +1 d ∂φ n + i,j,k ∆ t = 0 ⇒ φ n + i +1 ,j,k − φ n + i,j,k + φ n + i − ,j,k ∆ x + φ n + i,j +1 ,k − φ n + i,j,k + φ n + i,j − ,k ∆ y + φ n + i,j,k +1 − φ n + i,j,k + φ n + i,j,k − ∆ z + A n +1 xi + ,j,k − A nxi + ,j,k − A n +1 xi − ,j,k + A nxi − ,j,k ∆ t ∆ x + A n +1 yi,j + ,k − A nyi,j + ,k − A n +1 yi,j − ,k + A nyi,j − ,k ∆ t ∆ y + A n +1 zi,j,k + − A nzi,j,k + − A n +1 zi,j,k − + A nzi,j,k − ∆ t ∆ z = − e(cid:15) (cid:16) n ni + ,j + ,k + − n (cid:17) , (51) ∂S d ∂α n + i,j,k = ∂L n ∼ n +1 d ∂α n + i,j,k ∆ t = 0 ⇒ n n +1 i + ,j + ,k + + ∆ t ∆ x (cid:16) n n +1 i + ,j + ,k + v n +1 xi + ,j,k − n n +1 i − ,j + ,k + v n +1 xi − ,j,k (cid:17) + ∆ t ∆ y (cid:16) n n +1 i + ,j + ,k + v n +1 yi,j + ,k − n n +1 i + ,j − ,k + v n +1 yi,j − ,k (cid:17) + ∆ t ∆ z (cid:16) n n +1 i + ,j + ,k + v n +1 zi,j,k + − n n +1 i + ,j + ,k − v n +1 zi,j,k − (cid:17) = n ni + ,j + ,k + , (52) ∂S d ∂λ ni + ,j + ,k + = ∂L n ∼ n +1 d ∂λ ni + ,j + ,k + ∆ t = 0 ⇒ µ n + i,j,k = µ n − i,j,k − ∆ t ∆ x v nxi + ,j,k (cid:16) µ n − i +1 ,j,k − µ n − i,j,k (cid:17) − ∆ t ∆ y v nyi,j + ,k (cid:16) µ n − i,j +1 ,k − µ n − i,j,k (cid:17) − ∆ t ∆ z v nzi,j,k + (cid:16) µ n − i,j,k +1 − µ n − i,j,k (cid:17) , (53) ∂S d ∂µ n − i,j,k = ∂ ( L n − ∼ nd + L n ∼ n +1 d ) ∂µ n − i,j,k ∆ t = 0 ⇒ λ ni + ,j + ,k + + ∆ t ∆ x (cid:16) λ ni + ,j + ,k + v nxi + ,j,k − λ ni − ,j + ,k + v nxi − ,j,k (cid:17) + ∆ t ∆ y (cid:16) λ ni + ,j + ,k + v nyi,j + ,k − λ ni + ,j − ,k + v nyi,j − ,k (cid:17) + ∆ t ∆ z (cid:16) λ ni + ,j + ,k + v nzi,j,k + − λ ni + ,j + ,k − v nzi,j,k − (cid:17) = λ n − i + ,j + ,k + . (54)Eqs. (45)-(47) and (51) are constraints that restrict the solution manifold. Both initial-ization and evolution should be restricted by these numerical constraints. Eqs. (48)-(50)are explicit schemes used for updating the gauge field components. We emphasize that ina DEC framework and a rectangular lattice, the variational schemes (48)-(50) of Maxwell’s15quations are equal to the traditional FDTD method. Based on DEC, the electromagneticfields in a rectangular lattice are defined as, E n + xi + ,j,k = − A n +1 xi + ,j,k − A nxi + ,j,k ∆ t , (55) E n + yi,j + ,k = − A n +1 yi,j + ,k − A nyi,j + ,k ∆ t , (56) E n + zi,j,k + = − A n +1 zi,j,k + − A nzi,j,k + ∆ t , (57) µ H nxi,j + ,k + = A nzi,j +1 ,k + − A nzi,j,k + ∆ y − A nyi,j + ,k +1 − A nyi,j + ,k ∆ z , (58) µ H nyi + ,j,k + = A nxi + ,j,k +1 − A nxi + ,j,k ∆ z − A nzi +1 ,j,k + − A nzi,j,k + ∆ x , (59) µ H nzi + ,j + ,k = A nyi +1 ,j + ,k − A nyi,j + ,k ∆ x − A nxi + ,j +1 ,k − A nxi + ,j,k ∆ y , (60)where the temporal gauge φ n + i,j,k = 0 has been adopted explicitly. Then substitutingEqs. (55)-(60) into Eqs. (48)-(50), we obtain, E n + xi + ,j,k = E n − xi + ,j,k − ∆ t(cid:15) J nxi + ,j,k + ∆ t(cid:15) ∆ y (cid:16) H nzi + ,j + ,k − H nzi + ,j − ,k (cid:17) − ∆ t(cid:15) ∆ z (cid:16) H nyi + ,j,k + − H nyi + ,j,k − (cid:17) , (61) E n + yi,j + ,k = E n − yi,j + ,k − ∆ t(cid:15) J nyi,j + ,k + ∆ t(cid:15) ∆ z (cid:16) H nxi,j + ,k + − H nxi,j + ,k − (cid:17) − ∆ t(cid:15) ∆ x (cid:16) H nzi + ,j + ,k − H nzi − ,j + ,k (cid:17) , (62) E n + zi,j,k + = E n − zi,j,k + − ∆ t(cid:15) J nzi,j,k + + ∆ t(cid:15) ∆ x (cid:16) H nyi + ,j,k + − H nyi − ,j,k + (cid:17) − ∆ t(cid:15) ∆ y (cid:16) H nxi,j + ,k + − H nxi,j − ,k + (cid:17) , (63)where ( J x , J y , J z ) = ( env x , env y , env z ) is the current density. Eqs. (61)-(63) have the usualform of the standard FDTD method constructed in a Yee lattice [11]. It is well knownthat the standard FDTD is symplectic, which means the schemes (48)-(50) of Maxwell’s16quations are also symplectic. We will give a brief proof of this corollary in the following.Eq. (53) is an explicit scheme used for updating Lin’s constraint field. Eqs. (44), (52), and(54) are implicit schemes of the electron density and Lagrangian multipliers. It can be seenthat Eqs. (44)-(47), (52), and (54) make up an implicit cubic nonlinear algebraic systemthat can be solved for updating the electron density, velocity components, and multipliers.Effective and efficient nonlinear iterative methods and linear solvers are needed to solve thenonlinear algebraic equations. The complete iteration for the variational schemes is shownin Fig.3. FIG. 3. The complete iteration for the variational schemes of the hydrodynamic-electrodynamicmodel.
A very good property of the variational scheme is the conservation of Lagrangian structuregenerated by the discrete Lagrangian (43) L n ∼ n +1 d ( q n , q n +1 ), where q n = ( n nJ , v nJ , A nJ , φ n + J , α n + J , λ nJ , µ n − J )and the subscript J traverses all lattice points. We define two 1-forms, θ + n ∼ n +1 = ∂∂q n +1 i L n ∼ n +1 d (cid:0) q n , q n +1 (cid:1) · d q n +1 i , (64) θ − n ∼ n +1 = − ∂∂q ni L n ∼ n +1 d (cid:0) q n , q n +1 (cid:1) · d q ni , (65)which form a partition of the exterior derivative of the discrete Lagrangian d L n ∼ n +1 d = θ + n ∼ n +1 − θ − n ∼ n +1 . Then the Lagrangian noncanonical structure can be given as a closed2-form Ω n ∼ n +1 d = d θ + n ∼ n +1 = d θ − n ∼ n +1 as, Ω n ∼ n +1 d = ∂ ∂q ni ∂q n +1 j L n ∼ n +1 d (cid:0) q n , q n +1 (cid:1) · d q ni ∧ d q n +1 j . (66)17his discrete geometric structure is symplectic. By using the discrete dynamical equations(44)-(54) and taking the exterior derivative of action S d , we obtain, d S d = θ + N ∼ N +1 − θ − ∼ . (67)Eq. (67) means that one more exterior derivative of this 1-form leads to, Ω N ∼ N +1 d = Ω ∼ d . (68)The hydrodynamic-Maxwell system generated by Eq. (10) is naturally equipped with acontinuous Lagrangian symplectic structure. The conservation of discrete Lagrangian sym-plectic structure indicates the discrete variational principle generates a self-consistent finite-dimensional dynamical system that is a good analog of the continuous system. We shouldemphasize that preserving the discrete symplectic structure in the extended system (with La-grangian multipliers) is not a necessary and sufficient condition for preserving the geometricstructures in the original system. When the dynamical system is completely integrable, theKolmogorov-Arnold-Moser (KAM) theorem states that a weak perturbation in the Hamilto-nian will not break the invariant tori of the system, which is a theoretical root to ensure thegood numerical properties of the structure-preserving algorithms [21]. When it comes to ageneral dynamical system, although the KAM theorem is not always available, many numer-ical experiments show that the structure-preserving algorithms still exhibit good behaviorsin long-term simulations [18, 21].It can be directly examined that the gauge and translation symmetries are also preservedin the discrete dynamical equations. By introducing an arbitrary 0-form ψ , we can definethe discrete gauge transformation, φ (cid:48) n + i,j,k = φ n + i,j,k − ψ n +1 i,j,k − ψ ni,j,k ∆ t , (69) A (cid:48) nxi + ,j,k = A nxi + ,j,k + ψ ni +1 ,j,k − ψ ni,j,k ∆ x , (70) A (cid:48) nyi,j + ,k = A nyi,j + ,k + ψ ni,j +1 ,k − ψ ni,j,k ∆ y , (71) A (cid:48) nzi,j,k + = A nzi,j,k + + ψ ni,j,k +1 − ψ ni,j,k ∆ z . (72)Substituting Eqs. (69)-(72) into the discrete Lagrangian density (38), and summing over it ona universal discrete space-time manifold, we can obtain the gauge-invariant discrete action.Based on the discrete Noether’s theorem, it gives the discrete charge conservation law [30].18dditionally, the discrete Lagrangian density (38) is space-time-coordinate-independent,which means that the discrete action is translation-invariant. Based on the discrete Noether’stheorem, the Lagrangian momentum maps preserve the general energy-momentum density[30]. The advantages of variational schemes in long-term simulations are ensured by theconservation of the discrete symplectic 2-form, gauge, and translation symmetries, whichlead to the conserved quantities having long-term conservation and accuracy in simulations.The variational schemes constructed here are recognized as a structure-preserving Eu-lerian algorithm, which is convenient to implement and parallelize. There are also othertypes of structure-preserving algorithms constructed for hydrodynamic-electrodynamic sys-tems, such as the symplectic Smoothed-Particle-Hydrodynamics (SPH) method used forsimulating the double-fluid model of plasmas [51]. The symplectic SPH method is a Hybrid-Eulerian-Lagrangian (HEL) algorithm which avoids constraints and is suitable for simulatingplasma waves. But the metaparticle interpolation is tedious and time-consuming. The Eu-lerian form of the variational schemes in this work without tedious interpolation is moresuitable for nano-optics, although the constraints may have a singular performance undera few conditions. As the Lagrangian multipliers are pure mathematical variables, their ini-tialization is only determined by Eqs. (45) -(47). After the physical variables are initializedself-consistently, we can obtain the initial conditions of multipliers and Lin’s constraint fieldby solving Eqs. (45)-(47). C. Algebraic solvers
To implement the variational scheme, the solving procedure for nonlinear algebraic equa-tions (44)-(47), (52), and (54) is a core technique. In this work, we introduce the Newton-Krylov-type methods as primary algebraic solvers [52]. As the shell of a Newton-Krylov-type method, the Newton-Raphson algorithm is used as the basic nonlinear iterative methodwhich approximates and corrects the algebraic system in every numerical step.The nonlinear equations (44)-(47), (52), and (54) can rewritten in matrix form as, F (cid:0) X n +1 (cid:1) = 0 , (73)where F is a well-defined nonlinear vector function, X = ( α J , v xJ , v yJ , v zJ , n J , λ J ) T , and the19ubscript J traverses all lattice points. Then the Newton-Raphson iteration is given as, X ∗ = X − J − F ( X ) F ( X ) , (74)where X ∗ indicates new variables, and the Jacobian J F should be updated in every it-eration step. The detailed Jacobian elements can be found in App. A. The criterion || X ∗ − X || / || X || < ε is used to cutoff the iteration, where ε is a specified sufficientlysmall value.During every iteration step, we face a Jacobian inversion, which means a linear algebraicmatrix equation needs to be solved. Based on the Krylov subspace theory, there are manyefficient linear solvers. For example, the Generalized-Minimum-Residual (GMRES) method,the Incomplete-Cholesky-Conjugate-Gradient (ICCG) method, and the BICGSTAB method,have been constructed to solve a large sparse matrix equation [53–55]. In this work, theBICGSTAB method is introduced as the basic linear solver, because of its efficiency, stability,and parallel ability [55]. An alternative approach to solve the linear equations generatedin nonlinear iteration is the Jacobian-Free Newton-Krylov (JFNK) method, which replaces J F X ∗ with J F X ∗ ≈ [ J F ( X + ξ X ∗ ) − J F ( X )] /ξ , where ξ is a small perturbation [52, 56].A good feature of the JFNK method is that the Jacobian-vector product can be probedapproximately without forming and storing the Jacobian elements. In this work, we just usethe Newton-BICGSTAB iteration method to solve the nonlinear system, as the Jacobiancan be exactly derived conveniently (App. A). The JFNK approximation will be taken intoconsideration to accelerate the simulation in future work. IV. NUMERICAL EXPERIMENTSA. Bulk Plasmon
The first numerical experiment implemented in this work is the one-dimensional (1-D)bulk plasmon oscillation. This simulation is a numerical benchmark that is used to verify thevariational code. In this 1-D simulation, the metal is specified as silver, which means the theplasma frequency ω p = 1 . × Hz (for more details, see Sec.IV B). Then the spatial stepis chosen as 0 . c/ω p , and the temporal step is determined by the Courant-Friedrichs-Lewy(CFL) constraint, where CF L = 0 .
5. The numerical simulation domain is a 5000 1-D lattice,and the Message Passing Interface (MPI) is used as a parallel strategy. At the initial time,20 random perturbation of gauge field is introduced as the initial condition. After 10000steps simulation, the dispersion relation can be reconstructed by taking the Fast FourierTransform (FFT) of the gauge field. Fig. 4 shows the evolution of the gauge field in thesimulation, where the plasmon oscillation with frequency ω p can be recognized obviously.Fig. 5 plots the numerical dispersion relation versus the analytical one. The numerical result(contour plot) has a good consistency with the analytical dispersion relation (solid line) inthe truncation region k ∈ [0.00057,1.4] × m − . It shows that the accurate linear responsehas been reached over a wide range of the spectrum, which means the simulated plasmonicsystem is physically correct. FIG. 4. The evolution of the gauge field in the simulation of bulk plasmon oscillation.
To illustrate the good property of variational schemes in secular simulations, we plot theevolution of numerical error of total Hamiltonian, which is shown in Fig. 6 (c). After a longterm simulation, the numerical error is bounded by a very small value without coherentaccumulation. From the sub-Fig. 6 (a) and (b), we can find that the perturbation energyis exchanged cycle by cycle with frequency ω p between the electron and field componentsduring the oscillation, which drives the plasmonic motion.The result of first numerical experiment provide us with a basic verification of the numer-ical code used for simulating plasmonic phenomena. The advantage in conservation shownby numerical error is a footstone of secular simulation for a nonlinear system.21 IG. 5. Bulk plasmon dispersion relations: the numerical result (contour plot) versus the analyticalone (solid line). (a)(b)(c)
FIG. 6. Numerical error of the Hamiltonian in the simulation: (a) the Hamiltonian of losslosslessfree electron gas; (b) the Hamiltonian of self-consistent gauge field; (c) the total Hamiltonian error.The unit of Hamiltonian is Joule.
B. SPPs
A typical SPPs configuration is the infinite interface between metal and air. When itcomes to silver, the average electron number density n ≈ . × m − , and then theplasma frequency ω p = (cid:112) n e /(cid:15) m = 1 . × Hz. Based on the Drude model, themetallic permittivity can be given as (cid:15) ( ω ) = 1 − ω p /ω . The SPPs dispersion relation canbe analytically given as k x = k (cid:112) (cid:15)/ ( (cid:15) + 1), where k is the vacuum wave number, andthen the SPR frequency ω sp = 0 . × Hz. Fig. 7 shows the 2-D SPP configuration22sed in the simulations. The numerical simulation domain is a 200 ×
150 (metal 200 × x = ∆ z is 0 . × π/k x , the temporal step ∆ t is determined by the CFLconstraint CF L = 0 .
5, and the total simulation time is 4000 steps. To effectively implementthe simulations, the MPI is used to parallel the code.
FIG. 7. The 2-D SPP configuration used in the simulations.
Fig. 8 shows the dispersion relation at given wave lengths obtained by simulations. Bycomparing the numerical results (circle marks) with analytical dispersion curve, we findthat the simulations can reproduce characteristic dispersion relation of SPPs from nearlight cone to SPR regions. The resonance of wave means super resolution which is a coreproperty of near field optics. The inset of Fig. 8 plots the wave lengths and decay lengthsin the simulation region. It shows that the SPPs wave length (black line) and metallicdecay length (red line) are always less than the vacuum wave length (dashed line), butthe air decay length (blue line) is less than the vacuum wave length only in the near SPRregion ( ω > × Hz). The numerical results of SPPs wave lengths also have a goodprecision. Fig. 9 shows the slice of normalized SPPs field H y obtained by simulations, wherethe vacuum wave length λ = 240 nm. The Z-projection of Fig. 9 provide us with the SPPsmode structure on different sides of the silver-air interface. The localization of wave meanssubwavelength field structure and localized energy enhancement which can be used in many23 IG. 8. SPPs dispersion relations. Analytical vs numerical. Numerical results obtained at givenwave lengths have a good consistency with the analytical curve. The inset plots the wave lengthsand decay lengths at simulation region, which illustrates the sub-wavelength and energy localizationof SPPs.FIG. 9. The slice of normalized SPPs field H y obtained by simulations. The vacuum wave length λ = 240 nm. The insets plot the X- and Z- projections of H y at the region [35,65] × [25,75]. Itshows that the electromagnetic field of SPPs is strongly bounded by the interface, and the modehas a sub-wavelength spatial structure. branches of nano-electrodynamics. Just as the first numerical experiment, the conservedquantities, such as the total Hamiltonian, in these simulations are still bounded by thediscrete conservation laws.The numerical experiments implemented in this work verify the numerical schemes and24how the good properties of these schemes in secular simulations. It can be expected that thiswork will provide us with a powerful numerical tool for the first-principle based simulationstudy of plasmonics. V. CONCLUSION
In this work, we have constructed a class of variational schemes for the hydrodynamic-electrodynamic model of lossless free electron gas in quasi-neutral background to implementhigh-quality simulations of the SPPs. The Lagrangian density of lossless free electron gaswith a self-consistent electromagnetic field was established, and the conservation laws withconstraints were obtained. By using the DEC-based discretization, we reconstructed a fieldtheory on the discrete space-time manifold and constructed the variational schemes by dis-cretizing and minimizing the action. The variational schemes are nonlinear semi-explicit,which means the nonlinear solvers are needed. We introduced a hybrid Newton-BICGSTABmethod to solve the nonlinear algebraic equations involved in the variational schemes. In-stead of discretizing the partial differential equations, the variational schemes have superiornumerical properties in secular simulations, as they preserve the discrete Lagrangian sym-plectic structure, charge, and general energy-momentum density conservations. Two types ofnumerical experiments, i.e., bulk plasmon oscillation and 2-D SPPs, were implemented. Thenumerical results can reproduce characteristic dispersion relations of both types of plasmonicphenomena. The numerical errors of conserved quantities, e.g., the total Hamiltonian, inall examples are bounded by a small value after long-term simulations, which shows the ad-vantages of variational schemes in secular simulations. The variational schemes constructedfor the hydrodynamic-electrodynamic model can be used as a powerful numerical tool inplasmonics. Improvements, such as unstructured lattices, high-quality boundaries, more ef-ficient algebraic solvers, and loss and interband transition of real metal, will be developedin our future work.
ACKNOWLEDGMENTS
Q. Chen would like to thank H. Qin and J. Xiao for considerable help in differentialgeometry and field theory at University of Science and Technology of China. We also thank25he anonymous reviewer, whose criticism and suggestions strongly improved this paper. Thiswork is supported by the National Nature Science Foundation of China (NSFC-11805273,51477182), the CEMEE State Key Laboratory Foundation (CEMEE-2018Z0104B), and theTesting Technique Foundation for Young Scholars. Numerical simulations were implementedon the Shenma supercomputer at Institute of Plasma Physics, Chinese Academy of Sciencesand the TH-1A supercomputer at National Super Computer Center in Tianjin.
Appendix A: The Jacobian of Nonlinear Equations
To implement the Newton-Raphson iteration, we should calculate the Jacobian of nonlin-ear equations in every iteration step. Here we give the detailed Jacobian elements. Assumingthe nonlinear function f ( X n +1 ) is defined by Eq. (44), we obtain, ∂f ∂α n + i,j,k = 1 , (A1) ∂f ∂v n +1 xi + ,j,k = − m ∆ tv n +1 xi + ,j,k − e ∆ tA n +1 xi + ,j,k + ∆ t ∆ x (cid:16) α n + i +1 ,j,k − α n + i,j,k (cid:17) , (A2) ∂f ∂v n +1 yi,j + ,k = − m ∆ tv n +1 yi,j + ,k − e ∆ tA n +1 yi,j + ,k + ∆ t ∆ y (cid:16) α n + i,j +1 ,k − α n + i,j,k (cid:17) , (A3) ∂f ∂v n +1 zi,j,k + = − m ∆ tv n +1 zi,j,k + − e ∆ tA n +1 zi,j,k + + ∆ t ∆ z (cid:16) α n + i,j,k +1 − α n + i,j,k (cid:17) . (A4)Assuming the nonlinear function f ( X n +1 ) is defined by Eq. (45), we obtain, ∂f ∂v n +1 xi + ,j,k = mn n +1 i + ,j + ,k + , (A5) ∂f ∂n n +1 i + ,j + ,k + = mv n +1 xi + ,j,k + eA n +1 xi + ,j,k − α n + i +1 ,j,k − α n + i,j,k ∆ x , (A6) ∂f ∂λ n +1 i + ,j + ,k + = − µ n + i +1 ,j,k − µ n + i,j,k ∆ x . (A7)26ssuming the nonlinear function f ( X n +1 ) is defined by Eq. (46), we obtain, ∂f ∂v n +1 yi,j + ,k = mn n +1 i + ,j + ,k + , (A8) ∂f ∂n n +1 i + ,j + ,k + = mv n +1 yi,j + ,k + eA n +1 yi,j + ,k − α n + i,j +1 ,k − α n + i,j,k ∆ y , (A9) ∂f ∂λ n +1 i + ,j + ,k + = − µ n + i,j +1 ,k − µ n + i,j,k ∆ y . (A10)Assuming the nonlinear function f ( X n +1 ) is defined by Eq. (47), we obtain, ∂f ∂v n +1 zi,j,k + = mn n +1 i + ,j + ,k + , (A11) ∂f ∂n n +1 i + ,j + ,k + = mv n +1 zi,j,k + + eA n +1 zi,j,k + − α n + i,j,k +1 − α n + i,j,k ∆ z , (A12) ∂f ∂λ n +1 i + ,j + ,k + = − µ n + i,j,k +1 − µ n + i,j,k ∆ z . (A13)Assuming the nonlinear function f ( X n +1 ) is defined by Eq. (52), we obtain, ∂f ∂v n +1 xi + ,j,k = ∆ t ∆ x n n +1 i + ,j + ,k + , (A14) ∂f ∂v n +1 xi − ,j,k = − ∆ t ∆ x n n +1 i − ,j + ,k + , (A15) ∂f ∂v n +1 yi,j + ,k = ∆ t ∆ y n n +1 i + ,j + ,k + , (A16) ∂f ∂v n +1 yi,j − ,k = − ∆ t ∆ y n n +1 i + ,j − ,k + , (A17) ∂f ∂v n +1 zi,j,k + = ∆ t ∆ z n n +1 i + ,j + ,k + , (A18)27 f ∂v n +1 zi,j,k − = − ∆ t ∆ z n n +1 i + ,j + ,k − , (A19) ∂f ∂n n +1 i + ,j + ,k + = 1 + ∆ t ∆ x v n +1 xi + ,j,k + ∆ t ∆ y v n +1 yi,j + ,k + ∆ t ∆ z v n +1 zi,j,k + , (A20) ∂f ∂n n +1 i − ,j + ,k + = − ∆ t ∆ x v n +1 xi − ,j,k , (A21) ∂f ∂n n +1 i + ,j − ,k + = − ∆ t ∆ y v n +1 yi,j − ,k , (A22) ∂f ∂n n +1 i + ,j + ,k − = − ∆ t ∆ z v n +1 zi,j,k − . (A23)Assuming the nonlinear function f ( X n +1 ) is defined by Eq. (54), we obtain, ∂f ∂v n +1 xi + ,j,k = ∆ t ∆ x λ n +1 i + ,j + ,k + , (A24) ∂f ∂v n +1 xi − ,j,k = − ∆ t ∆ x λ n +1 i − ,j + ,k + , (A25) ∂f ∂v n +1 yi,j + ,k = ∆ t ∆ y λ n +1 i + ,j + ,k + , (A26) ∂f ∂v n +1 yi,j − ,k = − ∆ t ∆ y λ n +1 i + ,j − ,k + , (A27) ∂f ∂v n +1 zi,j,k + = ∆ t ∆ z λ n +1 i + ,j + ,k + , (A28) ∂f ∂v n +1 zi,j,k − = − ∆ t ∆ z λ n +1 i + ,j + ,k − , (A29) ∂f ∂λ n +1 i + ,j + ,k + = 1 + ∆ t ∆ x v n +1 xi + ,j,k + ∆ t ∆ y v n +1 yi,j + ,k + ∆ t ∆ z v n +1 zi,j,k + , (A30) ∂f ∂λ n +1 i − ,j + ,k + = − ∆ t ∆ x v n +1 xi − ,j,k , (A31)28 f ∂λ n +1 i + ,j − ,k + = − ∆ t ∆ y v n +1 yi,j − ,k , (A32) ∂f ∂λ n +1 i + ,j + ,k − = − ∆ t ∆ z v n +1 zi,j,k − . (A33)It can be seen that the Jacobian J F is a large sparse matrix in every iteration step. [1] R. Ritchie, Phys. Rev. , 874 (1957).[2] W. Barnes, A. Dereux, and T. Ebbesen, Nature , 6950:824 (2003).[3] E. Ozbay, Science , 5758:189 (2006).[4] R. Zia and M. L. Brongersma, nature nanotechnology , 426 (2007).[5] E. Hendry, F. Garcia-Vidal, L. Martin-Moreno, J. Rivas, M. Bonn, A. Hibbins, and M. Lock-year, Phys. Rev. Lett. , 123901 (2008).[6] S. Kim, J. H. Jin, Y. J. Kim, I. Y. Park, Y. Kim, and S. W. Kim, Nature , 757 (2008).[7] H. Atwater and A. Polman, Nature Materials , 205 (2010).[8] A. Boltasseva and H. A. Atwater, Science , 6015:290 (2011).[9] V. K. Valev, Langmuir , 15454 (2012).[10] M. Fang, Z. Huang, W. E. I. Sha, and X. Wu, IEEE J. Multiscale and Multiphys. Comput.Techn. , 194 (2017).[11] K. S. Yee, IEEE Trans. Antennas Propag. , 302 (1966).[12] R. F. Harrington, Field Computation by Moment Methods (MacMillan, New York, 1968).[13] A. Taflove and M. E. Brodwin, IEEE Trans. Microw. Theory Tech. , 623 (1975).[14] A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House Publisher, Boston, 1995).[15] Jr. J. D. Anderson,
Computational Fluid Dynamics: The Basics with Applications (McGraw-Hill, New York, 1995).[16] Q. Chen and B. Chen, Phys. Rev. E , 046704 (2012).[17] K. Feng, The Proceedings of 1984 Beijing Symposium on Differential Geometry and Differen-tial Equations (Science Press, Beijing, 1985) p. 42.[18] K. Feng and M. Qin,
Symplectic Geometric Algorithms for Hamiltonian Systems (Springer-Verlag, New York, 2010).
19] G. Benettin and A. Giorgilli, J. Statist. Phys. , 1117 (1994).[20] S. Reich, SIAM J. Numer. Anal. , 1549 (1999).[21] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration: Structure-PreservingAlgorithms for Ordinary Differential Equations (Springer, New York, 2002).[22] Y. K. Wu, E. Forest, and D. S. Robin, Phys. Rev. E , 046502 (2003).[23] E. Hairer, J. Sci. Comput. , 67 (2005).[24] S. A. Chin, Phys. Rev. E , 037701 (2009).[25] H. Qin, J. Liu, J. Xiao, R. Zhang, Y. He, Y. Wang, Y. Sun, J. W. Burby, L. Ellison, andY. Zhou, Nucl. Fusion , 014001 (2016).[26] M. Tao, Phys. Rev. E , 043303 (2016).[27] P. J. Morrison, Phys. Plasmas , 055502 (2017).[28] Q. Chen, H. Qin, J. Liu, J. Xiao, R. Zhang, Y. He, and Y. Wang, J. Comput. Phys. , 441(2017).[29] P. J. Morrison, Rev. Mod. Phys. , 467 (1998).[30] J. E. Marsden and M. West, Acta Numerica , 357 (2001).[31] A. Lew, J. E. Marsden, M. Ortiz, and M. West, Int. J. Numer. Meth. Engr. , 153 (2004).[32] M. West, Variational Integrators , Ph.D. thesis, California Institute of Technology (2004).[33] H. Qin and X. Guan, Phys. Rev. Lett. , 035006 (2008).[34] J. Li, H. Qin, Z. Pu, L. Xie, and S. Fu, Phys. Plasmas , 052902 (2011).[35] J. Squire, H. Qin, and W. M. Tang, Phys. Plasmas , 052501 (2012).[36] J. Squire, H. Qin, and W. M. Tang, Phys. Plasmas , 084501 (2012).[37] J. Xiao, J. Liu, H. Qin, and Z. Yu, Phys. Plasmas , 102517 (2013).[38] R. Zhang, J. Liu, Y. Tang, H. Qin, J. Xiao, and B. Zhu, Phys. Plasmas , 032504 (2014).[39] A. B. Stamm and B. A. Shadwick, IEEE Trans. Plasma Sci. , 1747 (2014).[40] B. A. Shadwick, A. B. Stamm, and E. G. Evstatiev, Phys. Plasmas , 055708 (2014).[41] J. Xiao, J. Liu, H. Qin, Z. Yu, and N. Xiang, Phys. Plasmas , 092305 (2015).[42] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendr¨ucker, J. Plasma Phys. , 905830401(2017).[43] Y. Shi, J. Xiao, H. Qin, and N. J. Fisch, Phys. Rev. E , 053206 (2018).[44] W. Newcomb, Nucl. Fusion Suppl. part2 , 451 (1962).[45] F. Sahraoui, G. Belmont, and L. Rezeau, Phys. Plasmas , 1325 (2003).
46] A. N. Hirani,
Discrete Exterior Calculus , Ph.D. thesis, California Institute of Technology(2003).[47] R. Hiptmair, Numer. Math. , 265 (2001).[48] D. N. Arnold, R. S. Falk, and R. Winther, Acta Numer. , 1 (2006).[49] D. N. Arnold, R. S. Falk, and R. Winther, Bull. New Ser., Am. Math. Soc. , 281 (2010).[50] M. Holst and A. Stern, Found. Comput. Math. , 263 (2012).[51] J. Xiao, H. Qin, P. J. Morrison, J. Liu, Z. Yu, R. Zhang, and Y. He, Phys. Plasmas ,112107 (2016).[52] P. N. Brown and Y. Saad, SIAM J. Sci. Stat. Comput. , 450 (1990).[53] Y. Saad and M. H. Schultz, SIAM J. Sci. Stat. Comput. , 856 (1986).[54] J. A. Meijerink and H. A. van der Vorst, Math. Comput. , 148 (1977).[55] H. A. van der Vorst, SIAM J. Sci. Stat. Comput. , 631 (1992).[56] D. A. Knoll and D. E. Keyes, J. Comput. Phys. , 357 (2004)., 357 (2004).