An Approximate Method for Modeling Plasmons on Graphene with Time- and Space-dependent Properties

1. Introduction

• The method targets current density for a propagating surface plasmon manipulated through time- and space-dependent Drude-weight changes.The authors motivate efficient computational methods as potential tools for designing photonic devices that use plasmonic phenomena on 2D materials.Source: §1 · pp. 1–2
• The introduction situates the work among graphene plasmonic device research, including integrative plasmonics, reconfigurable antennas, absorbers, and sensors.The cited examples include designs in and a proposed multi-band perfect absorber in.Source: §1 · p. 2
• The approach builds on and, which model time-dependent and then time- and space-dependent Drude weight, respectively.The present work continues.Source: §1 · p. 2
• Unlike the perturbational approach in, this work drops its assumption that mixing among different wave-number components is limited.That earlier assumption had received some numerical justification.Source: §1 · p. 2
• The model is framed as one-dimensional wave propagation with time- and space-dependent material properties, a setting the authors describe as sparsely analyzed in the open literature.The authors note that relates time-only and space-only 1D wave equations, allowing tools for propagation in inhomogeneous media to be applied.Source: §1 · p. 3

2. Model Equations

• The model uses Maxwell’s equations for a 2D transverse magnetic mode, with the graphene sheet lying along the x-axis and fields governed by jump conditions at the sheet.The electric permittivity and magnetic permeability are constant.Source: §2 · pp. 3–4
• On the sheet, the current density follows Drude’s law, with Drude weight D(x, t) and damping factor τ.The Drude weight varies in both space and time.Source: §2 · p. 4
• The initial fields and current correspond to the plasmon solution for constant Drude weight D0, with its wave number and decay parameter specified by the constant-weight case.The wave number ξ0 is a problem parameter; the time-frequency is obtained from a quartic equation and the decay parameter satisfies a dispersion relation.Source: §2 · pp. 4–5
• The current-density PIDE is used instead of solving the full field initial-value problem, with the stated motivation of avoiding 2D space truncation and absorbing boundary conditions.The PIDE formulation was derived in.Source: §2 · p. 5
• The Drude weight is modeled as a constant background plus a small perturbation, with α ≪ 1 and the perturbation bounded in the stated L2([0,T]), H1(R) norm.The paper pursues a perturbative approach under this assumption.Source: §2 · p. 5

3. Perturbation Expansion

• The Fourier-transformed current density is expanded as a series in α, with the constant-weight solution as its zeroth-order term.The Drude-weight transform is decomposed into a constant-background term and an α-scaled perturbation.Source: §3 · p. 6
• Substituting the expansion into the governing equation and grouping powers of α produces separate equations for each perturbation order.The paper treats the first-order equation and then extends the procedure to higher orders.Source: §3 · pp. 6–7
• The higher-order corrections are formulated as integral equations whose forcing is known from the preceding correction, enabling recursive computation.The integral equation is described as essentially a Volterra integral equation [19].Source: §3 · pp. 7–8

4. Solution method

• For each wave number ξ, the method solves the perturbation integral equations and inverts the Fourier transform to recover the spatial current corrections.The integral equations are solved over a finite time interval.Source: §4 · p. 8

4.1. Solvability of the integral equation

• If the perturbation is continuous in time, the integral equations have unique continuous solutions on [0, T] for each ξ.The solvability argument uses continuity of the kernel and forcing, with higher-order corrections handled inductively.Source: §4.1 · p. 8

4.2. Invertibility of the Fourier transform

• The higher-order invertibility argument proceeds inductively, with the second-order step requiring the perturbation in L2((0, t), H1(R)).Under that condition, the forcing term for the second-order equation is in L2((0, t), L2(R)).Source: §4.2 · pp. 11–12
• For perturbations in L2((0, t), H1(R)), every correction has a Fourier inverse for finite t and belongs to L2((0, t), L2(R)).For the first-order correction, the inverse exists for each finite t and has finite space-time L2 norm.Source: §4.2 · pp. 10, 12

5. Discretization

• For each Fourier mode, time discretization turns the integral equations into a matrix system, and inverse FFT recovers the spatial current density.The spatial interval is sampled uniformly and the inverse Fourier transform is approximated by an inverse DFT.Source: §5 · p. 13

5.1. Approximation Errors

• The kernel and integral operators are approximated with the trapezoidal rule, with stated errors of O(∆t^2).The kernel is approximated for l ≤ m, and the resulting operators are replaced by discrete counterparts.Source: §5.1 · pp. 13–14
• The linear-system solution has an l∞ error of O(∆t^2), while the total error analysis assumes no aliasing and includes O(∆t^2N) and O(∆ξ^2) terms.The analysis identifies three error terms; the final term requires a separate estimate.Source: §5.1 · p. 14

5.2. Computational Complexity

• The total computational cost is O(NM^3) + O(MN log N), with the O(NM^3) term likely to dominate in practice.The first term covers solving for the Fourier-domain values; Fourier inversion costs O(MN log N).Source: §5.2 · p. 15

6. Numerical Examples

• The MATLAB examples assume d(x, t) ∈ L2((0, T), H1(R)), use N = 1024 points on [−4π, 4π], and display [−π, π].The perturbation is chosen using smoothed rectangle terms with width a and sharpness r.Source: §6 · pp. 15–16

6.1. Propagating current density perturbation

• With ξ0 = ξ1 = 4 and ω0 = ω1, the correction terms are small, but the second-order contribution grows as time increases.The first- and second-order corrections are shown in Figure 1.Source: §6.1 · p. 16, Figure 1 · p. 16

6.2. Standing current density perturbations

• For ω0 = −ω1, standing waves arise in the correction terms, which remain small and concentrated within the support window of d(x, t).The second-order correction is included alongside the first-order term in the reported example.Source: §6.2 · p. 16, Figure 2 · p. 17

6.2.1. Growing current density perturbation

• Setting ξ0 = ξ1 = 4 and ω1 = √D0ξ0−ω0 produces a growing wave in the correction terms.Figure 3 shows the first- and second-order corrections for this modulation choice.Source: §6.2.1 · p. 17, Figure 3 · p. 17

7. Solution for Special Form of Drude Weight.

• For Drude weight D(x,t)=D0+αd(x,t), where d(·,t) is not in L2(R), the section presents a method to solve for the first-order perturbation approximation for j1(x,t).Source: §7 · p. 17
• For the special modulation, the transformed equation has right-hand-side support only at ξ=ν+ and ξ=ν−, and is solved analytically using the Laplace transform.The frequencies are ν+=ξ0+ξ1 and ν−=ξ0−ξ1.Source: §7 · p. 18

7.1. Inverting ˆ˜ja(s)

• Inverting ˜ja(s) requires contours for its poles and branch points; the apparent growing residue from p4 cancels with p2, leaving oscillatory or decaying residues, while the integrals are evaluated numerically.The branch points occur at s=±icν+.Source: §7.1 · pp. 19, 20, 22

7.2. Inverting ˆ˜jb(s)

• Inverting ˜jb(s) requires special cases because a pole may lie outside the contour.The formulas are provided in Appendix A.Source: §7.2 · p. 22

7.3. Inverting the Fourier Transform

• After Laplace inversion, the Fourier inversion is linear, so only the delta functions are inverted to recover the current terms.The resulting ˜ja(t) and ˜jb(t) do not depend on ξ.Source: §7.3 · p. 22

7.4.1. Traveling perturbational current density

• When ξ0=ξ1 and ω0=ω1, Figure 5 shows a traveling wave arising in the perturbation correction term; the solution should be compared with Figure 6 in.Source: §7.4.1 · p. 23

7.4.2. Standing perturbational current density

• With ξ0=ξ1 and ω0=−ω1, the perturbation term forms standing waves.The paper compares this case with Figure 7 in.Source: §7.4.2 · p. 24

7.4.3. Growing perturbational current density

• For ξ1=√D0ξ0−ω0, the perturbation term is a traveling wave that grows in amplitude.The paper compares this case with Figure 8 in.Source: §7.4.3 · p. 24

• The first method reformulates the governing PIDE as wave-number-parameterized Volterra integral equations and establishes existence, uniqueness, and spatial regularity of the perturbation.Source: §8 · p. 24
• Its discrete solver has complexity O(NM^3 + MN log N) and is reported more efficient than direct discretization with absorbing boundary conditions.The method is described as highly accurate for bounded spatiotemporal perturbations.Source: §8 · p. 24
• The second method uses contour integration and quadrature to invert transforms for a special form of Drude-weight perturbation.Numerical examples demonstrate the methods, which the authors present as design and predictive tools for integrated photonic devices leveraging plasmons on 2D materials.Source: §8 · p. 25

Appendix A. Inversion of ˆ˜jb(s)

• The inversion separates four cases according to whether ν− and ω0 − ω1 are zero, with the resulting pole and branch-cut structure treated case by case.Source: Appendix A · pp. 27, 29
• When ν− = 0 and ω0 − ω1 = 0, B− = 0 and the transformed expression is simplified before inversion.The appendix then uses partial fractions and a Bromwich contour with branch-cut contributions.Source: Appendix A · pp. 27–28
• For the case with ν− ≠ 0 and ω0 − ω1 ≠ 0, the analysis accounts for branch points at ±icν− and poles, using a notch when a pole lies below a branch point.The contour arcs vanish as R → ∞, and the lower-branch-cut contributions are evaluated around the notch.Source: Appendix A · pp. 29–30
• In another case, the branch points at ±icν− remain while the pole at −iω is absent, changing the inverse Laplace transform.Source: Appendix A · p. 32