Modeling of Subwavelength Gratings: Near-Field Behavior

Subwavelength gratings, with period shorter than the incident wavelength, have garnered significant attention in the fields of photonics, optoelectronics, and image sensor technology. In this research, we delve into the scattering characteristics of these gratings by employing the 2-dimensional point dipole approximation. Additionally, we propose a version of perturbation theory that relies on Fourier decomposition to obtain analytical expressions for the near-field behavior. We validate our models using numerical techniques such as boundary and finite element analysis. Notably, we explore how parameters like the grating period and slit width affect field enhancements. We demonstrate that our models produce qualitatively accurate results even for narrow slits.


I. INTRODUCTION
Subwavelength gratings (SWGs), characterized by periods shorter than the wavelength of the radiation they interact with, have garnered growing interest and found significant applications across various fields.Their versatility has been demonstrated in enhancing the focusing of light beams [1], contributing to advances in silicon photonics [2], enabling the development of optoelectronic devices exploiting surface plasmons [3], serving as mirrors in vertical-cavity surface-emitting lasers [4], and playing a pivotal role in image sensor technology [5,6].These structures provide compact and cost-effective solutions, even within the visible spectrum, offer exciting prospects for manipulating the polarization properties of emitted radiation [7], and enhancing the performance of InGaAs infrared polarization imaging [8].
The incorporation of double SWGs has shown promise in creating ultra-narrow resonances in the optical response of systems, with amplification factors of up to 10 5 [9].Moreover, SWGs are very useful in hybrid multiplexing technology to improve optical communication capacity by combining multi-wavelengths, multi-modes, as well as dual-polarizations [10].
The electrodynamics theory of diffraction relies on Maxwell's equations and the decomposition of their solutions into eigenfunctions.This theory includes the Rigorous Coupled-Wave Analysis (RCWA), which has primarily been developed for configurations of high symmetry.
Some examples of such configurations include relief profiles with square-wave, triangular, or saw-tooth shapes [11].Additionally, other notable cases involve periodic arrays of parallel infinite circular cylinders [12][13][14].In the case of relief profiles and ribbons, the eigenfunctions used are plane waves, whereas for cylindrical arrays, Bessel or Hankel functions serve as a basis.
There are various numerical methods available for analyzing this problem.The current cutting-edge techniques include Finite-Difference Time-Domain (FDTD) [15], Discrete Dipole Approximation (DDA) [16,17], Boundary Element Method (BEM) [18], and Finite Element Method (FEM) [19].However, each of these methods comes with inherent limitations and some numerical errors that makes it difficult to compare these approaches.
Consequently, there is a growing demand for new analytical solutions that can aid in optimizing photonic structures and validating numerical simulations.
Since the SWG period is much smaller than the radiation wavelength, a common simplification is to treat the grating as an effective continuous medium layer with averaged dielectric permittivity.This approach provides accurate results for reflection, transmission, and extinction when assuming that the effective medium is anisotropic [20].However, this simple averaging method fails to capture intricate interference phenomena like resonances at the onset of new diffraction modes (Rayleigh frequencies) [21,22] and near-field effects, like the sub-diffraction light focusing [23].The goal of this study is to use analytical treatment of SWGs with cylindrical elements.We apply the 2D point dipole approximation, similar to the 3D model of the interaction between dielectric spheres by Markel [24].Another approach used to get analytical expressions is the perturbation theory, valid for the dielectric permittivity close to unity |ε − 1| ≪ 1.It differs from the quantum-mechanical Born series because of the specific boundary conditions in the electrodynamics.A rather similar approach was applied earlier to the scattering problem for two parallel cylinders [25].
Section II introduces the Point Dipole Approximation (PDA), enabling the derivation of simplified analytical formulas.Section III describes the perturbation series based on the Fourier series developed for the near-field calculations.Section IV is devoted to numerical modeling by BEM and FEM to determine the near-field intensity card for electric and magnetic field, and comparing it to analytical formulas.Finally, Section V provides a summary of the conclusions.

A. Polarizability
We are examining an infinite array of cylindrical structures with a radius of a and a periodicity of L. An incident electromagnetic wave is coming from below along the y axis in its positive direction, with our coordinate system oriented such that z aligns with the axes of these parallel cylinders and x is perpendicular to them.We are condensing our analysis to a fundamental unit cell, as depicted in Figure 1 within the lattice focuses on the study of the near-field behavior.
This approximation is highly justified for a sparsely distributed lattice, especially when the radius of a cylinder is much smaller than the lattice period.As depicted in Figure 1(b), in the case of a sparse lattice, the field affecting each cylinder is approximately uniform.
The Point Dipole Approximation (PDA) can also be applied, to some extent, for a densely packed lattice when the gap dimension δ between neighboring cylinders is relatively small, i.e., when δ = L − 2a ≪ L. However, a noticeable error arises when assessing the interaction of a "dimer," which is a pair of adjacent cylinders, because the field near a polarized cylinder deviates from a simple dipole-like behavior.In such cases, it becomes necessary to consider additional terms in the multipole expansion.
The polarizability of a cylinder aligns well with the PDA and ultimately results in an effective lattice polarizability coefficient that differs from that of an individual cylinder.At normal incidence, the phase of incident field is the same for all cylinders.Consequently, the two-dimensional dipoles excited in them have the same phase.When evaluating the nearfield effects from these dipoles in the subwavelength regime, our calculations can be treated in the electrostatic limit, disregarding variations in the phase of the individual dipoles at the observation point.
The electric field vector of the incident wave lies along x-axis where e x is the unit orth.Total field E acting on each cylinder is the sum of external field E 0 and the field of all neighboring dipoles E ′ : The dipole moment per unit length d is proportional to the total field E: d = αE, where α is the polarizability of dielectric cylinder in a homogeneous external field [26]: Here ε is the dielectric constant of the cylinder material.
A two-dimensional dipole with number n directed along the x-axis creates field in the neighborhood of "central" cylinder with number 0. In the sub-wavelength (electrostatic) limit, the total field from all other cylinders at the center of cylinder n = 0 (x = 0) is given by the sum: Hence we obtain for the field of other dipoles and the dipole moment induced: Using formula (1), we get In this context, we can refer to the proportionality factor as the "grating polarizability" Eq. ( 2) highlights that the resonant field enhancement in the slits occurs in regions where the denominator is small.This phenomenon corresponds to the lattice dipole plasmon resonance.Specifically, such a resonance can be observed in a lattice composed of metallic cylinders in the optical frequency range where the real part of the dielectric constant, Re ε, is less than zero.
To achieve a strong resonance, it is essential to consider a "heterogeneous" lattice configuration with narrow slits between the cylinders.Unlike a dimer, which is not accurately described in PDA, the lattice resonance is a collective effect involving all the cylinders.For closely spaced scatterers, PDA still introduces some error, although it is not as pronounced as it is for a dimer configuration.

B. Near field
In the near-field approximation the field of n-th dipole at the observation point with Cartesian coordinates x, y is described by the expression where the dot inside parentheses denotes the scalar product, r n = e x x + e y y, vectors e x,y are unit orthants of the coordinate axes, x n = x − nL, r n = x 2 n + y 2 .The total field E x from all cylinders of the lattice, directed along the x-axis, is described by the infinite sum: For the y-component of the cylinder field E y at the observation point, respectively, we obtain: Note that both sums, give the periodic function of x.
Summing up the series, we get: The total electric field intensity of the scattered radiation is given by a simple formula: From Eq. ( 5), we can observe that the intensity of the scattered radiation exhibits periodic variations along the x-axis.In contrast, along the y-axis, the intensity diminishes symmetrically and quickly.
Figure 2 displays the distribution of scattered field intensity near a cylinder, in PDA.
One can see saddle points near the center of a gap.The positions of saddle points are x = (n + 1/2)L and y = 0, according to Eq. ( 5).It's important to note that this model results in an overestimation of the field intensity near the point dipole itself.Due to this overestimation, we have not included the regions around points x = nL, n = 0, ±1, . . ., and y = 0 in the figure.

III. PERTURBATION THEORY A. Fourier Series
PDA model provides a reasonably accurate description of the scattered electric field in the outer region, outside the cylinders.However, inside each cylinder, the model predicts a singularity, which is a point where the field becomes infinite.Additionally, this model does not account for the scattered magnetic field, as it effectively sets the rotor (curl) of the magnetic field to zero.Nevertheless, in reality, the scattered magnetic field is not zero.To address these limitations and obtain a more accurate representation of the electric field, we propose to consider the following equation: Derived from Maxwell's equations, it describes field E in a medium with magnetic vacuum permeability µ=1 and with space-variable dielectric permittivity ε(r).
The scattering of p-wave is commonly described using a single Helmholtz equation for the magnetic field.However, solving a periodic problem for a magnetic field using the Fourier method can be challenging due to the discontinuity of its derivative at the boundary caused by the imposed boundary conditions.To address this issue and effectively solve the problem, a pair of equations for a two-dimensional electric field is employed below.These equations ensure that the electric field is not only continuous but also maintains continuity with its derivative, thus facilitating a more tractable solution.
In the case of equidistant parallel cylinders, function ε(r) is expressed in terms of the theta-functions: For p-wave E z = 0, and components E x (x, y) and E y (x, y) remain in the equation (6).
Represent the components, the periodic functions of x, as Fourier series: At normal incidence of the external field this representation is natural.Substituting the sums, multiplying by exp(−2πipx/L), and integrating over an elementary cell, Fig. 1(a), from −L/2 to L/2, we obtain the following equations for the amplitudes: = e ik 0 y , thus the electric field at zero order has the specific physical meaning of the external field falling from bottom to top onto the grating.Equations for the first-order amplitudes U (1)

B. Solutions
At p ̸ = 0 we find a solution decreasing at infinity: At p = 0 there is no decreasing, then use different approach.First, we take any Green's function and describe the general solution The convenient Green's function to be the one that is zero when y < ỹ.
where J 0 , J 2 are the Bessel functions of the zeroth and second order, respectively [27].Within the first order B (1) 0 = i(ε − 1)k 0 πa 2 /2L, i.e. practically coincide with A (1) 0 .At last, the magnetic field can be found as a rotor of electric field vector: Estimating from Eq. ( 14), we can determine the small perturbation parameter by comparing the first and zeroth orders, given as: This observation suggests that the method remains valid for SWG, even when |ε − 1| ∼ 1, while k 0 a ≪ 1.The expansion is particularly effective for sparse gratings when a ≪ L.
The first-order correction is illustrated in Fig. 3, which includes 221 harmonics (with |m| ⩽ 110) chosen based on their ability to minimize phase variations between adjacent cylinders (δφ ≪ 1).In these plots, you can observe the maximum values of the scattered electric field intensity at the interfaces between the dielectric and free space, both on the right and left, see subfigure (a).There are saddle points at the geometric centers of the slits, existence of saddle points, and the appearance of a minimum strip at small values of y.
We repeat calculations for the same period L, but with a narrow slit δ = 10 nm.FEM yields the saddle points, too.The field enhancement factor increases by approximately 6 times in this case, as shown in Fig. 6.V. CONCLUSIONS Subwavelength gratings have a wide range of applications in photonics and optoelectronics.In this study, we analyze a periodic system of parallel cylinders.We use the 2-dimensional point-dipole approximation and Born series.We have derived analytical expressions for the electric and magnetic components in the near field.To validate our analytical model, we performed numerical simulations that prove its accuracy.Remarkably, we identified the emergence of saddle points at the geometric centers of the slits.Numerical simulations confirmed this observation.It is shown that the dipole approximation provides the correct qualitative description of the behavior in the slit, while the perturbation theory offers an approximate quantitative representation over the entire (x, y) plane.The findings from this research will enhance understanding facilitates optimized designs and pave the way for their application in cutting-edge optical technologies.

FIG. 1 :
FIG. 1: (a) The cross-section of an elementary cell in the periodic grating with cylindrical elements.(b) The periodical chain of point dipoles denoted by circles.

FIG. 2 :
FIG. 2: The electric field intensity distribution |E| 2 of the scattered radiation of the grating modeled by PDA (5) at L = 10 (in d 2 /L 4 units).In this representation, the white ovals highlight areas excluded from the diagram since their intensity exceeds 0.1.