Analyzing Near-Field Intensity Distribution
in Subwavelength Gratings
through Cylindrical Wave Decomposition

Abstract

The investigation into the scattering of plane waves by a periodic array of parallel cylinders utilizes the method of cylindrical wave decomposition, thereby reducing the problem complexity to a series of linear algebraic equations. This methodology proves particularly efficacious when the diameter of cylinders is significantly less than the wavelength of incident wave, resulting in a rapid diminution of the solution coefficients as a function of azimuth numbers. Such a reductionist approach facilitates the computation of scattered radiation intensity in near field. Subsequent cross-validation with numerical results corroborates the theoretical findings, showcasing a qualitative concordance between the two. This study underscores the efficacy of cylindrical wave decomposition in simplifying and accurately modeling wave scattering phenomena in structured media.

[label1]organization=Institute of Automation and Electrometry, Siberian Branch, Russian Academy of Sciences, addressline=1 Koptjug Avenue, city=Novosibirsk, postcode=630090, country=Russia [label2]organization=Novosibirsk State University, addressline=2 Pirogov Street, city=Novosibirsk, postcode=630090, country=Russia [label3]organization=Rzhanov Institute of Semiconductor Physics, Siberian Branch, Russian Academy of Sciences, addressline=13 Lavrentiev Avenue, city=Novosibirsk, postcode=630090, country=Russia

1 Introduction

Subwavelength gratings (SWGs), characterized by their periods being shorter than the wavelength of the interacting radiation, pose a significant challenge in the realm of nanoscale light-matter interactions [1]. Moreover, SWGs find extensive application across a broad spectrum of technological fields. They are integral to the functioning of optical filters, waveguides, and laser systems, facilitating the directed propagation of light [2]. Additionally, these gratings act as highly sensitive sensors capable of detecting chemical and biological substances by monitoring alterations in diffraction patterns upon analyte interaction [3], and they are instrumental in the operation of high-resolution spectral devices [4]. SWGs are also pivotal in the synthesis of metamaterials endowed with unique optical properties, including negative refraction indices and light superconductivity [5]. Their application extends to the field of lithography for the production of micro- and nanoelectronics, where they are utilized in the creation of light masks and patterns [6], and they play a significant role in the advancement of silicon photonics [7]. Furthermore, SWGs contribute to the enhancement of display technologies by elevating image quality and reducing power consumption [8], and they augment the efficiency of solar cells through increased light absorption, optimized polarization properties, and minimized reflection [9, 10]. In the domain of security and anti-counterfeiting, SWGs are deployed in the fabrication of protective holograms and optical elements that present challenges in replication [11]. The myriad applications of subwavelength diffraction gratings highlight their vast potential across diverse scientific and technological fields.

In far-field analysis, subwavelength gratings (SWGs) are modeled as a uniform layer with average thickness. In specific scenarios, the effective dielectric permittivity tensor must exhibit anisotropy [12]. An alternative analytical method proposes a planar grating model within an optically homogeneous environment, abstracting grating wires as a periodic dielectric constant variation [13]. Conversely, the near-field scenario presents intricacies, with detailed distributions. Scholarly discourse extensively explores cylinders, particularly the complexities of calculating the lattice sum [15, 14, 16, 17]. Recent studies of near-field distribution have unveiled phenomena like nanojets, indicating localized field intensity enhancements [18], and lower mode propagation within the grating [19]. Despite advancements, near-field distribution understanding during scattering remains limited. This aspect has been analyzed using various numerical methods, yet not incorporating cylindrical wave decomposition (CWD) [20]. Our research aims to elucidate intensity distribution using CWD and rigorous coupled-wave analysis, enhancing understanding of near-field dynamics in SWGs.

Section 2 introduces the cylindrical wave decomposition (CWD) method and concentrates on deriving algebraic equations for the coefficients of decomposition. When considering the SWG limit, we simplify the system to derive approximate formulas for these coefficients. Section 3 details the calculation of the magnetic field intensity using these derived formulas. Section 4 offers a comparative analysis of the results against those obtained from the Finite Element Method (FEM), perturbation theory, and the Discrete Dipole Approximation (DDA). Finally, Section 5 summarizes our findings.

2 Decomposition over cylindrical waves

Our study addresses the diffraction problem presented by a periodic array of parallel cylinders aligned along the $z$-axis, with the cross-sectional view illustrated in Fig. 1. We define the cylinder radius as $a$, its dielectric constant as $\varepsilon$, and the lattice period as $L$. We select a unit cell centered on the cylinder’s cross-section and investigate the scenario of normal incidence, with wave propagation along the positive $y$-axis.

Figure 1: Unit cell geometry: section in the plane perpendicular to the cylinder axis.

For the $p$-wave, we can write the Helmholtz equation for the magnetic field $H_{z}(x,y)$:

$$\bigtriangleup H_{z}+k^{2}H_{z}=0,$$ (1)

where $k=k_{0}\sqrt{\varepsilon}$ is the wave number inside the cylinder, $k_{0}=\omega/c$ is the wave number outside. Here $\omega$ is the frequency of radiation, $c$ is the speed of light in free space.

2.1 Boundary Conditions

The boundary conditions at the dielectric interface are simplified to the continuity of the tangential components of the electric and magnetic fields: $[E_{\tau}]=[H_{\tau}]=0$, with the square brackets indicating a discontinuity of the corresponding value inside. This results in two boundary conditions: the continuity of the magnetic field and its weighted normal derivative at the interface between the dielectric and free space at $r=a$:

$$H_{z}^{in}=H_{z}^{out},\quad\frac{1}{\varepsilon}\frac{\partial H_{z}^{in}}{\partial r}=\frac{\partial H_{z}^{out}}{\partial r},$$ (2)

where $r=\sqrt{x^{2}+y^{2}}$ is the radial coordinate, “in” indicates the interior region of the circle, and “out” refers to the exterior region. For the periodic structure problem, the periodicity condition on the cell’s sides complements Eq. (2):

$$H_{z}(-L/2,y)=H_{z}(L/2,y).$$ (3)

In cases of oblique incidence, the phase term $e^{ik_{0}xL}$ is introduced in condition (3) as per the Floquet theorem.

2.2 Fourier Series

The magnetic field in the observation point can be represented as the sum of incident and scattered wave $H_{z}^{s}$

$$H_{z}(\mathbf{r})=e^{ik_{0}y}+H_{z}^{s}(\mathbf{r}).$$ (4)

The incident plane wave in turn can be decomposed into a Fourier series in azimuth angle $\varphi$:

$$e^{ik_{0}y}=e^{ik_{0}r\sin\varphi}=\sum_{m=-\infty}^{\infty}e^{im\varphi}J_{m}(k_{0}r),\quad r>a.$$ (5)

The field inside $n$-th cylinder can also be written as an angular harmonic expansion:

$$H_{z,n}^{in}(\mathbf{r})=\sum_{m=-\infty}^{\infty}a_{m}J_{m}\left(k|\mathbf{r}-\mathbf{r}_{n}|\right)e^{im\varphi_{n}},\;|\mathbf{r}-\mathbf{r}_{n}|<a,\;\varphi_{n}=\arg(\mathbf{r}-\mathbf{r}_{n}).$$ (6)

Here $\mathbf{r}_{n}=(nL,0)$ are coordinates of the $n$th cylinder center. The scattered wave can be written as the sum of waves from each cylinder:

$$H_{z}^{s}(\mathbf{r})=\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}b_{m}H_{m}^{(1)}\left(k_{0}|\mathbf{r}-\mathbf{r}_{n}|\right)e^{im\varphi_{n}},\quad|\mathbf{r}-\mathbf{r}_{n}|>a.$$ (7)

Hereafter $\mathbf{r}=x+iy$ is 2-dimensional radius-vector. The coefficients of inner and outer series $a_{m},b_{m}$ are independent of the cylinder number, as follows from periodicity (3) and a specific excitation source in the form of a plane wave with normal incidence. The obtained Fourier series (7) is the required decomposition into cylindrical waves.

2.3 Equations for Coefficients

One can simplify the series exploiting the Graf addition theorem [21]

$$e^{i\nu\chi}Z_{\nu}^{(1)}(w)=\sum_{k=-\infty}^{\infty}e^{ik\alpha}Z_{\nu+k}^{(1)}(u)J_{k}(v),$$ (8)

where $Z_{\nu}$ is an arbitrary cylindrical function. Here, the vectors corresponding to the complex numbers $u,v,w$ form a triangle in the complex plane, with $v$ being the smallest of its sides, $\alpha,\chi$ are the interior angles of this triangle, opposite to sides $v,w$, respectively. It yields the Hankel function as a series of Bessel functions. Using Eq. (8) we get

$$e^{il\varphi_{n}}H_{l}^{(1)}(k_{0}|\mathbf{r}-\mathbf{r}_{n}|)=\sum\limits_{m=-\infty}^{\infty}e^{im\varphi}H_{(m-l)\sigma}^{(1)}(k_{0}r_{n})J_{m}(k_{0}r),$$ (9)

where $\sigma=\mbox{sign}\,n$.

Substituting (9) into (7) we find the scattered field

$$H_{z}^{s}(\mathbf{r})=\sum_{l=-\infty}^{\infty}b_{l}\left(\sum_{m=-\infty}^{\infty}J_{m}(k_{0}r)e^{im\varphi}F_{ml}+e^{il\varphi}H_{l}^{(1)}(k_{0}r)\right).$$ (10)

Here

$$F_{ml}=\sum_{n\neq 0}H_{(m-l)\sigma}^{(1)}(k_{0}r_{n})$$ (11)

is the lattice sum. Matrix element $F_{ml}$ is a superposition of cylindrical waves (Hankel functions) scattered by all the cylinders except of $n=0$:

$$F_{ml}=\left[\sum_{n=1}^{\infty}+\sum_{n=-\infty}^{-1}\right]H_{m-l}^{(1)}(k_{0}|n|L)=\sum_{n=1}^{\infty}H_{m-l}^{(1)}(k_{0}nL)\left(1+(-1)^{m-l}\right),$$ (12)

hence $F_{-m,-l}=F_{ml}$. The matrix elements $F_{ml}$ depend only on the difference of indices, i.e. the matrix is Teplitzian. It opens the possibility of fast algorithm of its inversion [22].

From the inner expansion (6), the outer expansion (10) and the boundary conditions (2), equating the coefficients at equal angular harmonics $m$, we find the algebraic system for the coefficients:

$$J_{m}+J_{m}\sum_{l}F_{ml}b_{l}+H_{m}b_{m}=a_{m}L_{m},$$ (13)

$$J_{m}^{\prime}+J_{m}^{\prime}\sum_{l}F_{ml}b_{l}+H_{m}^{\prime}b_{m}=\frac{a_{m}}{\sqrt{\varepsilon}}L_{m}^{\prime},$$ (14)

where $J_{m}\equiv J_{m}(k_{0}a),H_{m}\equiv H_{m}^{(1)}(k_{0}a),L_{m}\equiv J_{m}(ka)$.

We find the inner coefficient $a_{m}$ from (13), substitute them into (14), and get a system on the outer coefficients $b_{m}$:

$$b_{m}+C_{m}\sum_{l}F_{ml}b_{l}+C_{m}=0,$$ (15)

where

$$C_{m}=\frac{J_{m}L_{m}^{\prime}/\sqrt{\varepsilon}-J_{m}^{\prime}L_{m}}{H_{m}L_{m}^{\prime}/\sqrt{\varepsilon}-H_{m}^{\prime}L_{m}}$$ (16)

are auxiliary values, the coefficients of decomposition for a single isolated cylinder [23]. Azimuth number $m=\pm 1$ corresponds to dipole radiation; $m=0$ is responsible for isotropic scattering. The other terms of the expansion are small if the phase runup on the characteristic transverse dimension of the cylinder is small $ka\ll 1$. Fig. 2 confirms that the dipole scattering is the most significant for $p$-wave.

Figure 2: Coefficients $|C_{m}|$ at $\varepsilon=2.25$ and different $m=0,\dots,4$. From top to bottom $k_{0}a=0.3$, 0.2, 0.1.

3 Intensity Distribution

Here we find the near-field intensity distribution of the scattered wave with CWD. In the next section we compare our pattern with the results of calculation by other methods.

The scattered field in the external region is found from Eq (10):

$$H_{z}^{s}(\mathbf{r})=\sum_{n=-N_{\max}}^{N_{\max}}\sum_{m=-M_{\max}}^{M_{\max}}b_{m}e^{im\varphi}H_{m}^{(1)}(k_{0}|\mathbf{r}-\mathbf{r}_{n}|).$$ (17)

Coefficients $b_{m}$ are determined from the truncated system (15), with the summation constrained to $|m|\leqslant M_{\max}$ and $|n|\leqslant N_{\max}$. In our calculations, we set $M_{\max}=1$ and $N_{\max}=50$. The results indicate that the contribution of $m=2$ is two orders of magnitude smaller than that of $m=1$. To accommodate increased geometric filling factor $2a/L\leqslant 1$, a greater number of wires, $2N_{\max}+1$, must be included in the calculations.

Figure 3: Intensity distribution $|H_{z}(x,y)|^{2}$ at wavelength $\lambda=1.512~{}\mu$m, $N_{\max}=50,M_{\max}=1,\varepsilon=2.25,L=0.15~{}\mu$m: $a=0.01$ (a), $0.03$ (b), $0.05$ (c).

The distribution, as depicted in Fig. reff:CWD, is illustrated for a constant period $L=0.15mu$m across varying cylinder radii: $a=0.01,0.03,0.05mu$m, thereby modifying the filling factor. Distinct features of this distribution include the intensity $|Hz|^{2}$ peaking at the cylinder’s top and bottom edges, notably at the interface between free space and the dielectric at $x=nL,y=pma$. The peak’s amplitude escalates by an order of magnitude from one subfigure to the next, correlating with an increase in the filling factor. The peaks manifest as crescents, which shorten as the grating becomes sparser. Near $y=0$, a trough is observed, expanding at the midpoint of the gap between adjacent cylinders, with its width enlarging as the filling factor decreases. This widening gap is associated with a diminution in maximum intensity. The parameters $k0a=0.04,0.12,0.2$ suggest that for each cylinder configuration in Fig. reff:CWD(a), a sparse grating results in a pattern similar to that of a single isolated cylinder citeHarrington61. As the filling factor increases in Fig. reff:CWD(c), the interaction between cylinders becomes more pronounced.

4 Comparing with other methods

Figure 4: Intensity distribution $|H_{z}(x,y)|^{2}$ calculated at the same parameters, as in Fig. 3(a,b,c): finite elements method (a), perturbation theory (b), discrete dipole approximation (c).

4.1 Finite Elements Method

To compare the results, we used the COMSOL Multiphysics 6.1 environment [24] with the Wave Optics package. Geometry is defined as rectangle domain with wide $20\;\mu$m and height $1.8\;\mu$m containing a chain of cylinders. The domain material is air with $\varepsilon=1$, surrounded by layers of depth $d=0.3\;\mu$m where the perfectly matched layers boundary condition is specified. Chain is linear array of 100 identical cylinders (dielectric constant $\varepsilon=2.25$) displaced along $x$-axis with respect to parameters of problem: period $L=0.15\;\mu$m between the centers, radii are varying $a=0.01,0.03,0.05\;\mu$m. Background field defined as plane wave expansion of Gaussian beam is incident along $y$ axis: width $w_{0}=7.5\;\mu$m, focal plane at $y=0$, wavelength $\lambda=1.512\;\mu$m. We use inhomogeneous mesh of triangles ranging in size from 1 nm inside and in the immediate vicinity of the cylinder to 30 nm in the outer region.

4.2 Perturbation Theory

Due to the normal incidence of the external field on the periodic lattice of cylinders $\varepsilon(x,y)=\varepsilon(x+L,y)$, the scattered field will also be a periodic function $E^{s}(x,y)=E^{s}(x+L,y)$, which allows us to represent $E^{s}(x,y)$ as a Fourier series [20]:

$$E_{x}(x,y)=\sum_{n=-\infty}^{\infty}U_{n}(x)e^{2\pi inx/L},\quad E_{y}(x,y)=\sum_{n=-\infty}^{\infty}V_{n}(x)e^{2\pi inx/L}.$$ (18)

For $p$-wave we substitute series (18) into Maxwell’s equation for the electric field

$$\nabla\times\nabla\times\mbox{$E$}=\varepsilon(\mathbf{r})k_{0}^{2}\mbox{$E$}.$$ (19)

Equations for components are:

$$(E_{y})_{xy}-(E_{x})_{yy}=\varepsilon(\mathbf{r})k_{0}^{2}E_{x},\quad(E_{y})_{xx}+(E_{x})_{xy}=\varepsilon(\mathbf{r})k_{0}^{2}E_{y}.$$ (20)

Considering the right side of equations (19), (20) as a small perturbation, assuming the external field equal to $E_{x}^{(0)}=U_{0}^{(0)}=e^{ik_{0}y},E_{y}^{(0)}=0$, we obtain

$$V_{0}^{(1)}=0,\quad V_{n}^{(1)}=-\frac{2\pi in}{Lq_{n}^{2}}\frac{dU_{n}^{(1)}}{dy},\quad U_{0}^{(1)}=\frac{i(\varepsilon-1)k_{0}}{2L}\pi a^{2}e^{-ik_{0}y}$$

$$-(\varepsilon-1)k_{0}\int_{-a}^{a}\theta(y-\tilde{y})\sin\left(k_{0}(y-\tilde{y})\right)\frac{2\sqrt{a^{2}-\tilde{y}^{2}}}{L}e^{ik_{0}\tilde{y}}\,d\tilde{y},$$

$$U_{n}^{(1)}=-(\varepsilon-1)q_{n}^{2}\int_{-a}^{a}\frac{e^{-q_{n}|y-\tilde{y}|}}{2q_{n}}\frac{e^{ik_{0}\tilde{y}}}{\pi n}\sin\frac{2\pi n\sqrt{a^{2}-\tilde{y}^{2}}}{L}\,d\tilde{y}.$$ (21)

Here $q_{n}^{2}=(4\pi n)^{2}/L^{2}-k_{0}^{2}$ is the Fourier parameter of mode $n$, $V_{n}^{(1)}$ and $U_{n}^{(1)}$ are the principal terms of the perturbation ($V_{n}^{(m)}$ and $U_{n}^{(m)}$ are the $m$-th terms, respectively), $\theta$ is the Heaviside step function. If the scattered field being small, it is sufficient to retain $V_{n}^{(1)}$ and $U_{n}^{(1)}.$

Knowing $E_{x}^{((1))}$ and $E_{y}^{((1))}$, we can find magnetic field $H_{z}^{((1))}$:

$$H_{z}^{(1)}=\frac{1}{ik_{0}}\left[(E_{y})_{x}-(E_{x})_{y}\right]=\frac{1}{ik_{0}}\sum_{n\neq 0}k_{0}^{2}\frac{U_{n}^{(1)\prime}(y)}{q_{n}^{2}}-\frac{U_{0}^{(1)\prime}(y)}{ik_{0}}.$$ (22)

In contrast to expansion (18), series (22) converges faster due to the continuity of magnetic field in space (the electric field is piecewise continuous, discontinuities on the cylinder boundary worsen convergence). We hold $|n|\leqslant N_{\max}=30$ terms in Eq. (22), considering 61 terms of the series.

4.3 Discrete Dipole Approximation

The Discrete Dipole Approximation (DDA) method, extensively reviewed in [27] and originally introduced in [25, 26], has been a staple analytical tool in nanophotonics, aerosol and hydrosol optics, and astrophysics for decades. This method involves discretizing the scattering object into small cells treated as point dipoles, followed by solving the system of integral equations that describe dipole interactions. DDA is particularly adept at determining scattered fields and analyzing their dependence on incidence angle and wavelength.

While predominantly applied to three-dimensional problems, DDA has also shown efficacy in two-dimensional contexts, both in uniform [28, 29] and compound environments [30]. In this study, we employed 2-dimensional DDA to compute the scattered magnetic field from a finite chain of 49 cylinders excited by an incident plane wave. Each cylinder’s cross-section was divided into 352 square cells, equivalent to two-dimensional dipoles. By solving the equations for coupled dipoles, as detailed in [29], we obtained a total of 17,248 dipoles. These facilitated the calculation of scattering field in regions external to, and slightly distant from, the cylinder surfaces, due to the inaccuracy of dipole-emitted field calculations near edges.

4.4 Comparison

The intensity distribution derived from various computational methods — FEM, perturbation series, and DDA — is compared, utilizing the same parameters as for the CWD. As depicted in Fig. 4, analogous to Fig. 3, notable intensity peaks are observed at the cylinder tops and bottoms, alongside a relatively narrow trough along the $x$-axis at small $|y|$ values.

Fig. 4(a) illustrates results obtained through FEM within the COMSOL Multiphysics framework, revealing intensity peaks at the cylinder extremities and diminished value within the substantial gap between cylinders, akin to observations in Fig. 3(a). Fig. 4(b), derived using perturbation theory, exhibits an intensity increase correlating with the filling factor, mirroring the pattern seen in Fig. 3(b). Fig. 4(c) showcases the results from DDA, where white circles indicate uncalculated fields inside and immediately outside the cylinders, yet the external pattern closely aligns with the CWD results in Fig. 3(c).

Therefore, the characteristic features identified in CWD analyses are corroborated by findings from established numerical methods, affirming qualitative consistency across gratings with diverse filling factors.

5 Conclusions

The CWD method yields accurate near-field intensity distributions, consistent with computational electrodynamics and perturbation theory. Key characteristics of the distribution include: (i) Intensity peaks at coordinates $x=nL,y=\pm a$, where $n$ represents integers; (ii) A narrow trough near the $x$-axis, broadening with distance from between cylinders. In summary, CWD proficiently analyzes wave scattering in cylindrically symmetric systems, offering detailed near-field and modal insights.