Efficient method for the solution of Maxwell’s equations for nanostructured materials

The calculation of the electromagnetic field in nanostructured materials and nano-optoelectronic devices, when the wavelength of the incident radiation is comparable with the size of the structural elements, requires the exact solution of Maxwell's equations. In this case, a very promising numerical approach is the spectral element method, which combines the geometric flexibility of finite elements with high precision of spectral methods. In this paper the implementation of the spectral element method based on the Dirichlet-to-Neumann map for solving Maxwell’s equations is discussed. The application of the method for two-dimensional periodic structures, such as diffraction gratings and a metal nanowire array in a dielectric matrix, is demonstrated.


Introduction
The calculation of the electromagnetic field in nanostructured materials and nanooptoelectronic devices including image sensors, nanostructured solar cells, photonic crystals and diffraction gratings, when the wavelength of the incident radiation is comparable with the size of the structural elements, requires an exact numerical solution of Maxwell's equations. The spectral methods are one of the most powerful approaches for solving partial differential equations. According to these methods the solution of the problem is written as an expansion over known basis functions with subsequent calculation of expansion coefficients. When the permittivity -ε and the permeability -μ inside the domain change smoothly, global spectral methods provide exponentially fast convergence. That is, starting with the certain coefficient, the absolute values of the expansion coefficients decrease exponentially, and a few basis functions may be enough to achieve good accuracy. However, if the structure includes various materials the permittivity and permeability are not continuous functions and the convergence is only algebraic, in other words the absolute values of the coefficients decrease only as a certain power of the coefficient number. In this case, the spectral domain decomposition methods have the advantage, when the domain of interest is divided into smaller subdomains in each of which ε and μ are uniform and the decomposition of solution by local basis functions occurs separately in each of the subdomains. The subdomains are connected by using the interface conditions for electromagnetic fields. Thus, this approach combines geometric flexibility of finite elements with the high accuracy of spectral methods.
The first multi-domain spectral methods used the so-called strong form (patching) of the equations [1] that enforce the continuity of the solution and its first-order derivative at the interface between two adjacent subdomains. Later, in 1984, a spectral element method (SEM) was proposed that used the Galerkin weak formulation [2]. Now, some authors refer to spectral element methods only in this narrow sense. But due to apt turn of phrase another authors [3] use the term "spectral element methods" in wider sense for all methods that implement spectral approach inside the elements by analogy with finite element method. We also follow this interpretation here.
The main obstacle to the widespread use of SEM lies in much more complicated programming. To date, various techniques have been emerged to simplify the application of this method. In this paper we consider SEM based on the Dirichlet-to-Neumann (DtN) map. It can be used both to solve scattering problems and to search for eigenmodes. In the next part of the paper we shall demonstrate the application of the method to the two-dimensional periodic structures such as diffraction grating and metal nanowire array in a dielectric matrix; however, this approach can be extended to arbitrary three-dimensional calculations.

Method description
The general solution of Maxwell's equations for nonconical diffraction in the 2-D case can be represented as a linear combination of two fundamental polarizations: the transverse electric polarization (TE) and the transverse magnetic one (TM). In this paper, we devote attention to TM-case due to its higher complexity. Results for TE case are very similar.  , and  is the wavelength of the incident radiation. The relative permeability is set to be equal to unity. The incident wave is described by the ycomponent of the magnetic field grating. Note that using perfectly matched layers (PML) [4] this approach can be reformulated for a general aperiodic case. To implement spectral element method the domain of interest is divided by quadrilateral subdomains -spectral elements, so that the boundaries of different material (if any) coincide with the borders of the elements Fig. 1 [5]. If the subdomain boundaries are defined by curves , see Fig. 2, then the x variable will be written by new variables as: For the z variable instead of () If the length of one of the sides () i t Г goes to zero we obtain a curve-bounded triangle. In this case the transformation (2) also can be realized and the curve-bounded triangle will be mapped into the reference square. The derivatives are related by: By substituting (3) into (1) we get the equation to be solved for every subdomain. We shall connect the subdomains by utilizing a Dirichlet-to-Neumann (DtN) map. To calculate the matrix of the DtN map operator of a subdomain  , a rectangular grid of Following the influence matrix method [3,6], which is completely identical to the Schur complement method [6,7], a set of boundary conditions is solved sequentially, producing two square matrices NN  : in the first matrix F there will be N boundary conditions as columns, and in the second matrix  n F there will be columns of corresponding normal derivatives at the boundary. So, for an arbitrary boundary condition  f , the normal derivative   n f can be easily calculated: and taking into account boundary conditions on   , which for TM polarization are: ,, If the conditions on the boundaries 12 ,  are known, the conditions on the common boundary   will be: The strategy is similar to calculating values on a common boundary using the Smatrices [8]. Knowing the conditions on the boundary of the entire domain and by sequentially implementing Eq. (9) we can find the conditions on the boundaries of all subdomains.
If for two boundaries   f can be calculated: The described approach significantly simplifies the program implementation of spectral element method, allowing to divide the task into many small subtasks and to avoid operating with large matrices. The DtN map matrices of spectral elements can be calculated in parallel. It is also possible to connect the DtN maps of M individual elements into one in Log M stages, dividing the elements into groups, each consisting of two elements. In the Schur complement iterative method this problem is solved iteratively. Some iteration schemes are presented in [7]. The iterative algorithm is even more convenient for distributed parallel computing. In each computational node the DtN problem is solved for one or for a group of spectral elements and the nodes exchange only by the boundary condition information of the neighboring elements. In this point of view, the method is promising for calculations of large complex structures on distributed computing systems.
After the matrix of the DtN problem, that relates the values at the upper and lower boundaries of the domain (the entire structure) with the corresponding normal derivatives, is obtained, we can find the S-matrix, as well as the coefficients of the reflected and transmitted waves, depending on the coefficients of the incident wave. In front of and behind the domain assuming a medium where the dielectric constant is uniform or varies only along the x axis, the magnetic field can be written through the eigenfunctions () s x  of the one-dimensional problem, as is done in modal methods [9]: In particular, in a homogeneous medium with a dielectric constant l  : and S-matrix: one can obtain equation for S-matrix: If instead of the collocation method we use the Galerkin method with linearly independent test functions multiplying by corresponding weights, which are represented as rows of the matrices 12 ,  ΨΨ , the equation will be: In this case, the number of points on the boundary can be greater than the number of coefficients j c  , and it is convenient to use the zeros of the Legendre polynomial as grid points for more precise integration. For the case presented in Fig. 1  c waves from (16) and using (14), (12) and (9) the boundary conditions for all subdomains can be found. After that, by solving the Dirichlet problem in each subdomain, the magnetic field is calculated.
The spectral element method can be also used for iterative calculation of eigenmodes, for example, in resonators, or as an integral part of the modal methods. In paper [9] such iterative approach was used to calculate hundreds of thousands of eigenmodes within modal method based on 1-D spectral elements.

Numerical results
The first example of application of the method discussed in this paper is the grating sketched in Fig. 1a)   The wavelength λ of the incident TM polarized field forming an angle of incidence θ=30 0 with respect to the normal to the device plane is equal to 3 / 4 p . The squared magnetic field for the considered grating is showed up in Fig. 1b). Fig. 3 illustrates the error of the reflectance of the discussed grating as function of the total number of nodes N x along the axis x which is equal to the number of elements multiplying by the number of nodes in one element along this axis. For simplicity, we define the same number of nodes on all sides of each element. As reference value we use reflected zero order R 0 ≈0.15407111704. When we fix the number of elements and increase the number of nodes in every element (the line with red squares) the convergence is exponentially fast. That is the so-called "prefinement" strategy. When the number of nodes in every spectral element is fixed (green circles and blue diamonds lines) convergence is only algebraic and much slower. The second example is the array of metal nanowires (silver) with permittivity The permittivities of the external environment is ε 1 =1, the wavelength λ of the TM polarized incident field (with θ=30 0 ) is equal to 1.5 μm, / 3 / 2 p   . The radiuses of nanowires are 1 =45,