# Numerical analysis of discrete geometric method on plasmonic structures

Shuai Yan, Xiaoyu Xu, Christoph Pflaum, Zhuoxiang Ren

► To cite this version:

Shuai Yan, Xiaoyu Xu, Christoph Pflaum, Zhuoxiang Ren. Numerical analysis of discrete geometric method on plasmonic structures. IEEE Transactions on Magnetics, 2016, 52 (3), pp.7206304. 10.1109/TMAG.2015.2490278 . hal-01315771

HAL Id: hal-01315771

https://hal.sorbonne-universite.fr/hal-01315771

Submitted on 29 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical Analysis of Discrete Geometric Method on Plasmonic Structures

Shuai Yan<fnref n="1" />, Xiaoyu Xu<fnref n="1" />, Christoph Pflaum<fnref n="2" />, and Zhuoxiang Ren<fnref n="1" />,<fnref n="3" />

Discrete geometric method (DGM) is implemented for solving Maxwell's equations on plasmonic structures. Surface plasmons introduce non-derivable field components at the metallic/dielectric interface that can influence the accuracy of DGM. An analysis shows that the proper setting of material parameters at the interface can increase the accuracy of numerical solutions. Numerical examples with a structured mesh and an unstructured mesh provide further evidence for the fact that a proper interface treatment is important to obtain accurate results with relatively less degrees of freedom for the simulation of plasmonic effects with DGM. Coupling with its flexibility and explicit formulation, DGM can be used as an accurate and fast solver for problems involving complex plasmonic structures.

**Index Terms**—Computational electromagnetics, convergence of numerical methods, discrete geometric method (DGM), surface plasmons (SPs).

## I. INTRODUCTION

PLASMONICS are defined by the interacting processes between the electromagnetic waves and the conducting electrons at metallic interfaces or nanostructures, which leads to an enhanced optical near field of a sub-wavelength dimension. In recent years, accompanied with the development of nanofabrication techniques, a surface plasmon (SP) shows its great potential in a large variety of applications [1], including sub-wavelength optics, data storage, photocatalysis, solar cells, and biosensors.

Plasmonics can be explained by the classical Maxwell's framework. However, the accompanied evanescent fields form non-classical solutions of Maxwell's equations. In the other word, the parallel components of the electromagnetic fields are not differentiable at the metallic interface. To improve the accuracy of the numerical methods in modeling the evanescent fields, there are two common strategies. The first is to introduce analytical interface conditions in the discrete models [2], [3]. This can bring additional complexity to the formulation of the numerical solvers, especially for complex geometries. The second strategy is to introduce finer meshes in the vicinity of the interfaces. The field solver stays unchanged in this way, but the number of degrees of freedom (DoFs) increases. To save the computational resource in this case, before refining the discretization mesh with a brute force, it is also important to analyze how the discretization error arises at the interface, accordingly choose more proper ways for field discretization and material modeling at the interface, so that the requisite accuracy can be achieved with a relative smaller number of DoFs.

Discrete geometric method (DGM) is based on a differential geometry interpretation of Maxwell's system. Within this framework, different types of finite methods, including a finite-difference (FD) time-domain method, finite integrate techniques (FIT), and a finite-element method (FEM), can be unified in a general form [4], [5]

$$C^T h = d_t d + j \quad (1)$$

$$C e = -d_t b \quad (2)$$

$$d = M_{\varepsilon} e, b = M_{\mu} h, j = M_{\sigma} e \quad (3)$$

where e, h, d, b, and j are the line or surface integrals of electric field E, magnetic field H, electric flux density D, magnetic flux density B, and electric current density J, respectively, C is the discrete version of the curl operator, and d_t denotes the time derivative. With this algebraic form, (1) and (2) keep the metric-free nature of Faraday's law and Ampere's law. The metric information is condensed in the constitutive matrices M_ε, M_μ, and M_σ in (3). The three constitutive matrices actually give discretization to the continuous Hodge operator that maps 1-differential forms to 2-differential forms in E<fnref n="3" /> space. Thus, no matter how a discrete method is derived, it is only distinct from other methods in the manner to discretize the Hodge operator. DGM construct the constitutive matrices by a physically natural way [5], [6]. With diagonal constitutive matrices [6], the electromagnetic fields can be solved explicitly in the time or frequency domain [7], and therefore, the need to solve large-scale linear systems can be avoided. The method can be applied on an unstructured mesh, so that curved geometric features can be modeled more accurately.

In this paper, we derive the discretization error of DGM on a planar interface sustaining SPs. Based on the analysis, we propose a manner for specific field discretization and constitutive matrices construction to achieve better accuracy on a fixed mesh. Numerical examples are given

<footnote n="1">¹Institute of Microelectronics, Chinese Academy of Sciences, Beijing 100029, China</footnote>

<footnote n="2">²Department of Computer Science, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen 91058, Germany</footnote>

<footnote n="3">³Sorbonne Universités, University Pierre et Marie Curie, Paris 75005, France</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Demonstration of the time-harmonic SPP solution at a single metallic and dielectric interface. Red curve: distribution of the field component at the cross section marked by the black dashed line.</figcaption></figure>

for calculations on both the structured and unstructured meshes.

## II. SURFACE PLASMONS

There are two existing states of SPs: 1) SP polaritons (SPPs) are propagating evanescent waves along the metal/dielectric interfaces and 2) localized SP resonances (LSPRs) are non-propagating resonances of electrons at the metallic nanostructures.

The simplest geometry sustaining SPs is a single, flat interface between a dielectric, non-absorbing half space with a positive real dielectric constant $ε_d$ and an adjacent conducting half space with a dielectric function $ε_m(ω) = ε_{m,r}(ω) + iε_{m,i}(ω)$. The requirement of metallic character implies $ε_{m,r}(ω) < 0$ [1]. SPP state only exists for a transverse magnetic mode, and in this case, the analytical solution on this geometry can be obtained by solving the following Maxwell's equation of magnetic field for both half planes:

$$ \frac{1}{\varepsilon} \left( \frac{\partial^2 \tilde{H}_y}{\partial z^2} + \frac{\partial^2 \tilde{H}_y}{\partial x^2} \right) + k_0^2 \tilde{H}_y = 0. \quad (4) $$

The solutions are explicitly express as follows [1]:

$$ \left\{ \begin{aligned} \tilde{H}_y(x, y, z) &= A \exp(-k_d(z - z_0)) \exp(i\beta x) \\ \tilde{E}_x(x, y, z) &= -iA \frac{k_d}{\omega\varepsilon_0\varepsilon_d} \exp(-k_d(z - z_0)) \exp(i\beta x) \\ \tilde{E}_x(x, y, z) &= A \frac{\beta}{\omega\varepsilon_0\varepsilon_d} \exp(-k_d(z - z_0)) \exp(i\beta x) \end{aligned} \right. \quad z > z_0 \quad (5) $$

$$ \left\{ \begin{aligned} \tilde{H}_y(x, y, z) &= A \exp(k_m(z - z_0)) \exp(i\beta x) \\ \tilde{E}_x(x, y, z) &= iA \frac{k_m}{\omega\varepsilon_0\varepsilon_m} \exp(k_m(z - z_0)) \exp(i\beta x) \\ \tilde{E}_x(x, y, z) &= A \frac{\beta}{\omega\varepsilon_0\varepsilon_m} \exp(k_m(z - z_0)) \exp(i\beta x) \end{aligned} \right. \quad z \le z_0. \quad (6) $$

The tilde indicates that the fields are time-harmonic components of the original electric and magnetic fields. A is an arbitrary constant. The complex wave number $β$ of SPP takes the value $β = k_0(ε_mε_d/(ε_m + ε_d))^{1/2}$. $k_m$ and $k_d$ are related to the dielectric constants by an interface condition $k_m/k_d = -ε_m/ε_d$.

The distribution of field around the interface is shown in Fig. 1. The field is enhanced at the dielectric/metallic interface and decreases evanescently in both materials. For LSPRs, the field distribution possesses the same property at the surface of the metallic nanoparticles with the same interface conditions.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Demonstration of two types of field distribution at the dielectric/metallic interface. Case 1: electric fields lying on the interface. Case 2: cutting across the interface.</figcaption></figure>

## III. DISCRETE GEOMETRY METHOD AND DISCRETIZATION ERROR ESTIMATION

DGM is derived directly from the topology of the discrete mesh, without going through the process of discretizing a continuous differential or integral equation as with FD, FIT, or FEM. For DGM, the constitutive matrices are defined directly based on the geometry of the primal and dual meshes, where the electric field and the magnetic field are defined, respectively. Taking $M_ε$ as an example, it is a diagonal matrix, and each of its nonzero elements corresponds to a DoF e locating either on the edge $L_m$ of primal or dual cells. For an edge-surface dual pair ($L_m, A_m$) in homogeneous regions, the corresponding diagonal element of $M_ε$ is defined as $M_ε = ε|A_m|/|L_m|$. When either $L_m$ or $A_m$ cuts across two elements, $(M_ε)_{m \cdot m}$ is determined according to an equivalent circuit perspective [6]. Fig. 2 shows the two cases with $L_m$ lying on the material interface and cut by the material interface on a 2-D structured mesh. We denote the dual pairs as $(L_m^{(p)}, A_m^{(d)})$ and $(L_m^{(d)}, A_m^{(p)})$, respectively, and for the two cases, e is defined on the primal cell in the first case and on the dual cell for the second case.

The corresponding diagonal element of $M_ε$ is then calculated by

$$ (M_{\varepsilon})_{m \cdot m} = \frac{\varepsilon_d |A_1^{(d)}|}{|L^{(p)}|} + \frac{\varepsilon_m |A_2^{(d)}|}{|L^{(p)}|}, \quad (\text{Case 1}) \quad (7) $$

$$ (M_{\varepsilon})_{m \cdot m} = \frac{\varepsilon_d \varepsilon_m |A^{(p)}|}{\varepsilon_d |L_2^{(d)}| + \varepsilon_m |L_1^{(d)}|}, \quad (\text{Case 2}). \quad (8) $$

The definitions can be easily generalized to an unstructured mesh.

When the constitutive relation (3) is well defined for DGM, the question is about how the discrete system approximates the continuous Maxwell system arises, especially with the non-differentiable field components caused by SPs. In the following part, we evaluate the truncation error of the discrete equations when approximating to the continuous ones. Let $e^*$, $h^*$, $d^*$, $b^*$, and $j^*$ be the exact solutions of Maxwell's system. On the discrete mesh, they fulfill a linear system as follows:

$$ C^T h^* = d_t d^* + j^* \quad (9) $$

$$ C e^* = -d_t b^* \quad (10) $$

$$ d^* = M_{\varepsilon} e^* + R_{\varepsilon}, \quad b^* = M_{\mu} h^* + R_{\mu}, \quad j^* = M_{\sigma} h^* + R_{\sigma} \quad (11) $$
<table><caption>TABLE I
LOCAL TRUNCATION ERROR AT THE INTERFACE</caption><thead><tr><th></th><th>First term in (15)</th><th>Second term in (15)</th><th>Summation</th></tr></thead><tbody><tr><td>Case 1</td><td>O(h<sub>x</sub><sup>2</sup>) + O(h<sub>z</sub><sup>2</sup>)</td><td>O(h<sub>x</sub><sup>2</sup>) + O(1)</td><td>O(h<sub>x</sub><sup>2</sup>) + O(1)</td></tr><tr><td>Case 2</td><td>O(h<sub>x</sub><sup>2</sup>) + O(h<sub>z</sub>)</td><td>O(h<sub>x</sub><sup>2</sup>) + O(h<sub>z</sub>)</td><td>O(h<sub>x</sub><sup>2</sup>) + O(h<sub>z</sub>)</td></tr></tbody></table>

where $R_ε$, $R_μ$, and $R_σ$ are discretization errors on constitutive relations.

Let $\mathbf{r}^e = \mathbf{e}^* - \mathbf{e}$ and $\mathbf{r}^h = \mathbf{h}^* - \mathbf{h}$, and subtract (1)–(3) from (9)–(11). We obtain the following system:

$$C^T \mathbf{r}^h = i\omega M_\varepsilon \mathbf{r}^e + M_\sigma \mathbf{r}^e + i\omega \mathbf{R}_\varepsilon + \mathbf{R}_\sigma \quad (12)$$

$$C \mathbf{r}^e = -i\omega M_\mu \mathbf{r}^h + M_\sigma \mathbf{r}^e + i\omega \mathbf{R}_\mu. \quad (13)$$

Note that the time dependence is removed by the time-harmonic assumption. The variables take the same expressions without ambiguity. Solving out $\mathbf{r}^e$ from (12) and substituting it into (13) yield

$$\mathbf{r}^h - \omega^2 M_\mu^{-1} C \left( M_\varepsilon + \frac{M_\sigma}{i\omega} \right)^{-1} C^T \mathbf{r}^h = \mathbf{R}. \quad (14)$$

For the single interface SPP geometry explained in Section II, (14) is the equation of discretization error corresponding to (4), and the truncation error of the discrete equation to the continuous one is defined as

$$\mathbf{R} = M_\mu^{-1} \mathbf{R}_\mu - \omega^2 M_\mu^{-1} C \left( M_\varepsilon + \frac{M_\sigma}{i\omega} \right)^{-1} (i\omega \mathbf{R}_\varepsilon + \mathbf{R}_\sigma). \quad (15)$$

Assuming the calculation is carried on a 2-D computational domain with a single dielectric/metallic interface, and the domain is discretized by a uniform mesh with a mesh size $(h_x, h_z)$. We can evaluate $\mathbf{R}$ for each DoF $\mathbf{h}$ straightforwardly, since the constitutive matrices are diagonal. The error terms $\mathbf{R}_\varepsilon$, $\mathbf{R}_\sigma$, and $\mathbf{R}_\mu$ are calculated by the fact that the exact solutions [given in (5) and (6)] are smooth piecewise and fulfill the interface conditions at the interface. The exact solutions are expanded in Taylor series, and the low-order terms are canceled out to get the order of the error terms.

Given any DoF $\mathbf{h}$ locating at least one cell away from the interface, the local truncation error $\mathbf{R}$ can be easily evaluated as $O(h_x^2 + h_z^2)$. When the interface involves, we calculate the two subtraction terms of $\mathbf{R}$ in (15) for both cases shown in Fig. 2. The results are presented in Table I. For both cases, the error with respect to $h_z$ is dominant, since the interface is perpendicular to the z-axis. In case 1, the truncation error is a constant at the interface. Therefore, the discrete formulation is not consistent to the continuous differential equation in the $L_\infty$ sense. While in case 2, the consistency is fulfilled, and the convergence rate is of $O(h_z)$. In a weaker $L_2$ sense, a convergence with an order of $O(h_z^{1/2})$ can also be reached for case 1, since the local errors in the whole computational domain is averaged when evaluating the $L_2$ norm.

# IV. NUMERICAL EXAMPLES

## A. Geometry With a Single Planar Interface

For a further validation of the results from the theoretical analysis, we use DGM to solve an SPP problem at a planar

<figure><img src="image_4.png" /><figcaption>Fig. 3. Changes of $L_\infty$ and $L_2$ errors as the mesh is refined in the two cases of field distribution. $p$ is the corresponding convergent rate.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Geometry of the periodic Ag nanogratings.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Transmission spectra calculated with different mesh sizes for both cases. (a) Case 1 and (b) case 2 of field distribution.</figcaption></figure>

Ag/air interface. Both cases in Fig. 2 with electric components on the primal (case 1) and dual edges (case 2) are considered. The length in the x-direction of the computational domain is set as two periods of the traveling SPP. The height is 100 nm with an interface at $z_0 = 50$ nm with Ag below and air above. Imposing the Dirichlet boundary condition at the boundary of the computational domain, we calculate the numerical solution for different levels of discretization. The $L_\infty$ and $L_2$ errors of the numerical solution to the exact solution are shown in Fig. 3. The values are normalized according to their maximum. The convergent rate is defined as $p = \log_2 (||\mathbf{r}_N^h||/||\mathbf{r}_N^h||)$, where $N$ is the number of grid points. The results agree well with the theoretical analysis in Section III.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Magnetic field strength at a nano-Ag cylinder (unit: A/nm).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Absorption cross sections. (a) Results for a coarser mesh with 932 elements. (b) Results for a finer mesh with 2180 elements.</figcaption></figure>

## B. Periodic Nanogratings

We then consider a periodic silver (Ag) nanograting structure with geometry shown in Fig. 4. This kind of structures is of interest as transparent electrodes in optoelectronic devices because of their extraordinary transmittance, which is closely related to their LSPR properties [7], [8]. The computational domain is discretized with the uniform rectangular meshes. We set periodic boundary conditions for the sides and the perfect-matched layers (PMLs) [9] for the top and bottom of the domain. The calculation is also carried on for both cases with electric components lying on and cutting across the interface. The transmission of the structure is then calculated from the calculated electromagnetic fields. A reference curve from [8], which is calculated from the FEM is shown in Fig. 5(a) and (b) for a comparison. Results show that for coarse meshes with a mesh size of 5 and 2.5 nm, larger deviation can be observed in case 1 than in case 2. As the mesh is refined, the resulted spectra in both cases converge. This agrees with the theoretical analysis applied on the single metallic/dielectric interface. The results state again that with electric field defined on a dual mesh (case 2), we can achieve more accurate results with a coarser mesh.

## C. Single Nanocylinder

For curved geometries, DGM is well suitable for the simulation on an unstructured mesh. The dual mesh is constructed by connecting the circumcenters of the primal mesh [6]. We consider geometry with an infinite long Ag cylinder. The diameter of the Ag cylinder is 30 nm. The cylinder is placed in the center of a 120 nm × 120 nm domain. PMLs are set at all boundaries to truncate the computational domain. Fig. 6 shows the strength of the computed magnetic field for an incident wave with 340 nm wavelength on an unstructured mesh.

The calculation is carried out with interface settings for both cases shown in Fig. 2 on a coarser mesh with 932 elements and a finer mesh with 2180 elements. We calculate the absorption cross section of the cylinder as a quantity of interest. Results are shown in Fig. 7. The analytical values from the Mie calculation are given as a reference [10]. For both meshes, local setting in case 2 results in a better fit to the reference curve as expected.

## V. CONCLUSION

In this paper, we apply DGM on the simulation of plasmonic effects at metallic/dielectric interfaces. The theoretical and numerical analyses suggest that with electric fields cut across the metallic/dielectric interface, more accurate solutions can be obtained with a relatively coarser mesh when solving with the DGM formulations. Combining with its other advantageous properties, including geometric flexibility and explicit formulation, DGM can be further used as an accurate and fast numerical solver for an electromagnetic simulation on complex plasmonic structures.

## ACKNOWLEDGMENT

This work was supported in part by the National Science Foundation of China under Grant 51407181 and in part by the Beijing Natural Science Foundation under Grant 4144090.

## REFERENCES

[1] S. A. Maier, *Plasmonics: Fundamentals and Applications*. New York, NY, USA: Springer-Verlag, 2007, pp. 5–30.

[2] A. Mohammadi, T. Jalali, and M. Agio, "Dispersive contour-path algorithm for the two-dimensional finite-difference time-domain method," *Opt. Exp.*, vol. 16, no. 10, pp. 7397–7406, 2008.

[3] C. Classen, E. Gjonaj, U. Römer, R. Schuhmann, and T. Weiland, "Modeling of field singularities at dielectric edges using grid based method," *Adv. Radio Sci.*, vol. 9, no. 3, pp. 39–44, 2011.

[4] M. Clemens and T. Weiland, "Discrete electromagnetism with the finite integration technique," *Prog. Electromagn. Res.*, vol. 32, pp. 65–87, 2001.

[5] A. Bossavit, "Generalized finite differences' in computational electromagnetics," *Prog. Electromagn. Res.*, vol. 32, pp. 45–64, 2001.

[6] Z. Ren and X. Xu, "Dual discrete geometric methods in terms of scalar potential on unstructured mesh in electrostatics," *IEEE Trans. Magn.*, vol. 50, no. 2, Feb. 2014, Art. ID 7000704.

[7] S. Yan, J. Krantz, K. Forberich, C. Pflaum, and C. J. Brabec, "Numerical simulation of light propagation in silver nanowire films using time-harmonic inverse iterative method," *J. Appl. Phys.*, vol. 113, no. 15, p. 154303, 2013.

[8] J. Y. Lee, S. T. Connor, Y. Cui, and P. Peumans, "Solution-processed metal nanowire mesh transparent electrodes," *Nano Lett.*, vol. 8, no. 2, pp. 689–692, 2008.

[9] S. D. Gedney, "An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices," *IEEE Trans. Antennas Propag.*, vol. 44, no. 12, pp. 1630–1639, Dec. 1996.

[10] C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles*. New York, NY, USA: Wiley, 1983, pp. 57–81.