Symmetry Reduction in FEM Optics Modeling of Single and Periodic Nanostructures

Abstract

Numerical optics modeling is an invaluable tool in the design of nanostructures for nanophotonics applications where diffraction effects often lead to complex dependency between the nanostructure geometry and its optical properties and response. In order to analyze, design, and optimize such nanostructures, computationally efficient numerical optics modeling methods are required. One way to improve the numerical performance is to exploit symmetries found in many optics problems. By identifying equivalencies and restrictions arising from symmetry, it can be possible to simplify the problem at hand, which is the essence of symmetry reduction. However, applying symmetry reduction in optics modeling problems is not trivial. To the best of our knowledge, symmetry reduction has so-far been applied in finite element method (FEM) optics models only in those specific cases where an incident plane wave shares symmetries with the nanostructure geometry. In this work, we show how to extend the symmetry reduction of FEM optics models to the case of nonsymmetric plane-wave incidence, demonstrate such reduction with numerical examples of incident plane wave absorption in a single nanowire and a periodic nanowire array, and discuss the achieved gains in computational efficiency.

1. Introduction

Numerical optics modeling [1] is an invaluable tool in the design of nanostructures with the desired optical properties and response, e.g., scattering, absorption, and emission, for nanophotonics applications including solar cells [2,3], photodetectors [4,5], and LEDs [6]. Such nanostructures typically have feature sizes comparable to the considered wavelength of light and diffraction effects often lead to complex dependency between the nanostructure geometry and its optical properties and response [1,2,7,8,9]. Such problems can rarely be treated analytically, and experimenting with fabricated device prototypes is costly and time-consuming. Instead, numerical optics modeling can be used to predict these properties for the given geometry and materials and to analyze the underlying optical response further, e.g., by performing modal analysis or inspecting the spatial distribution of the fields inside the nanostructures. In the design optimization of nanophotonic structures for applications, optimization algorithms tend to require many iterations where the numerical models are solved. Overall, there is a clear need for computationally as efficient optics modeling methods as possible.

One way to enhance the numerical performance with the most popular numerical modeling methods, including the finite element method (FEM), the finite-difference time-domain (FDTD) method, and the Fourier modal method (FMM), is to exploit symmetries found in many optics problems. For example, nanostructures and nanostructure arrays can have geometries that effectively remain unchanged under certain mirroring and rotation operations, and the optical response, for given incident fields or for emission from within the system, can be assumed identical for symmetrical directions. Furthermore, symmetry can also be used for assessing which optical modes can couple to each other. By identifying equivalencies and restrictions arising from symmetry, it can be possible to simplify the problem at hand. This is the essence of symmetry reduction which is a well-known concept and has been used very successfully in other fields like computational chemistry (see e.g., ref. [10]).

However, applying symmetry reduction in optics modeling problems is not trivial. For example, in a problem with a field incident on a periodic array of nanostructures, the symmetries of the system depend on the array geometry, the geometry of nanostructures in the unit cell, and the configuration of the incident field. Inspecting and treating the symmetry properties of systems in a formal way can be achieved with the well-known tools found in group theory (see e.g., ref. [10]). Indeed, the properties of Maxwell’s equations under symmetry transformations [11,12] as well as exploiting point group symmetries in solving scattering problems, e.g., with boundary integral formulations [13,14], multiple-scattering formalism [15], discrete dipole approximation [16], and FMM [17], has been considered in the literature. However, applying symmetry reduction to similar optics problems in real-space-based FEM or FDTD implementations has so-far been lacking. We focus specifically on FEM since it is a versatile and popular method to obtain solution to light scattering problems for both single nanostructures and periodic nanostructure arrays (note however that the formalism presented below should be directly applicable for FDTD). To the best of our knowledge, symmetry reduction has so-far been applied in FEM optics models only in those specific cases where an incident plane wave shares symmetries with the nanostructure geometry.

In this work, we show how to extend the symmetry reduction of FEM models to work with nonsymmetric plane-wave incidences. We start with decomposing the field solution in the full simulation domain (the original problem, OP) into symmetry respecting modes, found using tools from group theory, as has been done also in the symmetry reduction of FMM optics modeling [17]. We then translate the symmetry modes to different boundary conditions in the smaller size, symmetry reduced simulation domain, resulting in several computationally less demanding subproblems (SPs). Finally, since we consider the linear Maxwell equations, we can obtain the OP solution by properly adding together the SP solutions. To demonstrate our symmetry reduction method, we provide numerical examples of incident plane wave absorption in a single nanowire and a periodic nanowire array and discuss the achieved gains in computational efficiency compared to models without symmetry reduction.

2. Methods

In our symmetry reduction method for FEM optics modeling, light is considered as an electromagnetic field governed by the Maxwell’s equations, and we focus on the case of scattering and absorption of incident light. More specifically, we assume time harmonic fields ($\mathbf{E}(\mathbf{r},t)=\Re \left\{\mathbf{E}\left(\mathbf{r}\right)e^{i\omega t}\right\}$, where $\omega =2\pi c_0/{\lambda }_0$ is the angular frequency, $c_0$ is the speed of light in vacuum, and ${\lambda }_0$ is the wavelength in vacuum), local and linear response, and media that are nonmagnetic (vacuum permeability ${\mu }_0$) and isotropic, such that they can be described by a spatially-varying, frequency-dispersive, complex-valued, scalar refractive index $n(\mathbf{r},\omega )=\sqrt{{\epsilon }_{\mathrm{r}}(\mathbf{r},\omega )}$, where ${\epsilon }_{\mathrm{r}}(\mathbf{r},\omega )$ is the relative permittivity. Furthermore, we consider source-free regions, i.e., fields far away from charge or current sources that give rise to the incident light, meaning that the incident light is considered as a boundary condition. With these assumptions, the Maxwell’s equations yield the well-known wave equation for the complex-valued electric field:

$$\nabla \times \nabla \times \mathbf{E}\left(\mathbf{r}\right)-k^2(\mathbf{r},\omega )\mathbf{E}\left(\mathbf{r}\right)=0,$$

(1)

where $k(\mathbf{r},\omega )=n(\mathbf{r},\omega )k_0=n(\mathbf{r},\omega )2\pi /{\lambda }_0$ is the wave number. The complex-valued magnetic field is then connected to the electric field by Maxwell’s equations as $\mathbf{H}\left(\mathbf{r}\right)=-1/\left(i\omega {\mu }_0\right)\nabla \times \mathbf{E}\left(\mathbf{r}\right)$. For the following discussion we note that Equation (1) is linear, i.e., any linear combination of solutions $\mathbf{E}\left(\mathbf{r}\right)$ is also a solution, and that the electric field can be solved independently from the magnetic field. We omit explicitly writing the frequency dependence for the rest of our discussion to keep the equations more concise.

2.1. Symmetries in Single Nanostructures and Periodic Arrays

In this work, we consider nanostructures on top of a semi-infinite substrate. Therefore, regardless of the nanostructure geometry, the system can have at most a single axis of rotational symmetry that is normal to the substrate surface. We choose the z-axis of our coordinate system to be parallel to this single principal axis of rotational symmetry (as is the usual convention). Note that, in the absence of the nanostructures, the system with a single interface between the semi-infinite superstrate and substrate media is continuously rotational symmetric with respect to this axis (the rotational symmetry axis is of order infinity). We also define the substrate–superstrate interface to coincide with the $z=0$ plane. Furthermore, the considered geometries can obviously only have vertical mirror symmetry planes, that is, planes with normal perpendicular to the z-axis, and cannot have the inversion symmetry (except for the specific case of a homogeneous surrounding where the substrate material is the same as the material on top of the substrate).

In the case of a single nanostructure, the symmetries can be determined by considering the $xy$ cross-section (note that the nanostructure geometry does not need to be translation invariant along the z-axis as long as the cross-section maintains its symmetries for varying z). In our numerical examples, we consider nanowire structures with regular hexagon cross-section as illustrated in Figure 1. The system consisting of the nanowire and the substrate hence belongs to the ${\mathrm{C}}_6\mathrm{v}$ point group with the following symmetry elements [10]: the six-fold rotational axis ${\mathrm{C}}_6$, three mirror planes coinciding with the hexagon vertices (${\sigma }_{\mathrm{v}}^1$, ${\sigma }_{\mathrm{v}}^2$, and ${\sigma }_{\mathrm{v}}^3$), and three mirror planes coinciding with the hexagon edge midpoints (${\sigma }_{\mathrm{d}}^1$, ${\sigma }_{\mathrm{d}}^2$, and ${\sigma }_{\mathrm{d}}^3$). These symmetry elements are illustrated in Figure 1b. The point group can be expressed with the corresponding symmetry operations as ${\mathrm{C}}_{6\mathrm{v}}=\{\widehat{\mathrm{E}},{\widehat{\mathrm{C}}}_6^1,{\widehat{\mathrm{C}}}_3^1,{\widehat{\mathrm{C}}}_2,{\widehat{\mathrm{C}}}_3^2,{\widehat{\mathrm{C}}}_6^5,{\widehat{\sigma }}_{\mathrm{v}}^1,{\widehat{\sigma }}_{\mathrm{v}}^2,{\widehat{\sigma }}_{\mathrm{v}}^3,{\widehat{\sigma }}_{\mathrm{d}}^1,{\widehat{\sigma }}_{\mathrm{d}}^2,{\widehat{\sigma }}_{\mathrm{d}}^3\}$, where $\widehat{\mathrm{E}}$ denotes the identity operation, ${\widehat{\mathrm{C}}}_n^m$ denotes the operation of rotation by m times the angle $2\pi /n$ around the rotational symmetry axis, and ${\widehat{\sigma }}_{\mathrm{v}}^1$, ${\widehat{\sigma }}_{\mathrm{v}}^2$, ${\widehat{\sigma }}_{\mathrm{v}}^3$, ${\widehat{\sigma }}_{\mathrm{d}}^1$, ${\widehat{\sigma }}_{\mathrm{d}}^2$ and ${\widehat{\sigma }}_{\mathrm{d}}^3$ denote the operations of mirroring across the mirror symmetry planes. We also define the x- and y-axis such that ${\widehat{\sigma }}_{\mathrm{v}}^1={\widehat{\sigma }}_{xz}$ and ${\widehat{\sigma }}_{\mathrm{d}}^2={\widehat{\sigma }}_{yz}$, meaning the operation of mirroring around the $xz$- and $yz$-plane, respectively.

Other examples of high-symmetry cross-sections include regular polygons and a circle. The nanostructures with regular polygon cross-sections belong to the point groups ${\mathrm{C}}_{N\mathrm{v}}$, where $N\ge 3$ is the number of vertices in the polygon. These point groups have the N-fold rotational symmetry axis and similarly defined vertical mirror symmetry planes as the hexagonal cross-section. For our symmetry reduction in the FEM modeling, it is important to note that all these point groups have the ${\mathrm{C}}_{\mathrm{s}}=\{\widehat{\mathrm{E}},{\widehat{\sigma }}_{xz}\}$ as a subgroup and, with even N, also the ${\mathrm{C}}_{2\mathrm{v}}=\{\widehat{\mathrm{E}},{\widehat{\mathrm{C}}}_2,{\widehat{\sigma }}_{xz},{\widehat{\sigma }}_{yz}\}$ group. Since ${\widehat{\mathrm{C}}}_2={\widehat{\sigma }}_{xz}{\widehat{\sigma }}_{yz}={\widehat{\sigma }}_{yz}{\widehat{\sigma }}_{xz}$, these groups are defined by the $xz$ and $yz$ mirror symmetry planes. Furthermore, there exist a multitude of other cross-sections that also have the ${\mathrm{C}}_{\mathrm{s}}$ or ${\mathrm{C}}_{2\mathrm{v}}$ either as their point group or as a subgroup thereof. For example, a rectangular or elliptical cross-section would have the ${\mathrm{C}}_{2\mathrm{v}}$ point group. The circular cross-section, on the other hand, is a special case as it has continuous rotational symmetry (rotational symmetry axis of order infinity). It is certainly possible to freely rotate the Cartesian coordinate system around the z-axis (the rotational symmetry axis) and consider the ${\sigma }_{xz}$ and ${\sigma }_{yz}$ mirror symmetries. However, in order to take full advantage of the rotational invariance, one should rather be able to symmetry reduce to the $xz$-plane with $x\ge 0$ and expand the fields in cylindrical harmonics. However, such symmetry reduction is not the focus of our work.

In the case of periodic arrays of nanostructures, symmetries of the geometry are defined by both the single nanostructure geometry and the 2D lattice to which they are arranged on the substrate surface. However, if the nanostructures have circular cross-section, the symmetries are defined entirely by the lattice arrangement. There are five distinct 2D Bravais lattices: monoclinic, tetragonal (square), hexagonal, orthorhombic (rectangular), and orthorhombic centered (rectangular centered). These lattices and their unit cells are illustrated in Figure 2. With the monoclinic lattice, we can only have the ${\mathrm{C}}_2$ symmetry element, but with all the other lattices, we can find primitive or non-primitive rectangular unit cells with both $xz$ and $yz$ mirror symmetry plane, provided that the nanostructure geometry also has these symmetries. Single nanostructure point groups and the five Bravais lattices combined result in a total of 17 distinct plane groups (see ref. [17] and the references therein for details). It is often possible to find rectangular unit cells in the resulting geometries that have the $xz$ or $yz$ mirror symmetry plane or both.

2.2. Well-Known Symmetry Reduction Methods with FEM

We find that in FEM optics modeling, exploiting symmetry to reduce the simulation domain size and hence the computational cost is generally applied only to a limited extent. In periodic structures with discrete translational symmetry, the structure can be reduced to a single unit cell with periodic boundary conditions. Typically, Floquet-type periodic boundary conditions are used when the unit cell has pairs of opposing boundaries and the field solution on these boundaries is known, due to the periodicity, to differ only by a phase factor, as given by the Floquet vector ${\mathbf{k}}_{\mathrm{F}}$. In the case of an incident plane wave, the Floquet vector is given by the in-plane components of the wave vector ${\mathbf{k}}_{\mathrm{i}}=k_{\mathrm{i}x}{\mathbf{u}}_x+k_{\mathrm{i}y}{\mathbf{u}}_y+k_{\mathrm{i}z}{\mathbf{u}}_z$ of the incident plane wave: ${\mathbf{k}}_{\mathrm{F}}=k_{\mathrm{i}x}{\mathbf{u}}_x+k_{\mathrm{i}y}{\mathbf{u}}_y$, where ${\mathbf{u}}_x$, ${\mathbf{u}}_y$, and ${\mathbf{u}}_z$ are the unit vectors in x-, y-, and z-direction. This type of symmetry reduction is routinely applied to most problems with periodicity (with the exception of, e.g., periodic structures that extend over only a few periods). However, exploiting mirror symmetries for a single nanostructure, or for the unit cell, is less often utilized.

A mirror symmetry plane corresponds to either a perfect magnetic conductor (PMC) or perfect electric conductor (PEC) condition when the electric field of the light is symmetric or antisymmetric with respect to the plane, respectively. The mirror symmetry operation corresponding to the plane will leave the electric field vector components parallel to the plane unchanged while reversing the sign of the component perpendicular to the plane. Therefore, if the electric field is symmetric with respect to the plane, then on the plane it must be parallel to the plane and, following from Maxwell’s equations, the magnetic field must be perpendicular to the plane. Hence, the PMC condition is fulfilled. Similarly, if the electric field is antisymmetric with respect to the plane, then on the plane it must be perpendicular to the plane with the magnetic field parallel to the plane, thus fulfilling the PEC condition. It is common in optics modeling based on Maxwell’s equations to first solve only for the electric field in the system and then obtain the magnetic field from that solution. Therefore, it is sufficient to consider the symmetries only for the electric field.

To the best of our knowledge, mirror symmetry reduction in FEM modeling has so-far been restricted to using PMC or PEC boundaries when the electric field (either incident or emitted) fulfills these conditions, i.e., when the field shares these mirror symmetries of the geometry. There are two common situations when this applies. First, when a plane wave is incident on a single nanostructure such that the plane of incidence coincides with a mirror symmetry plane of the geometry, the simulation domain can be cut in half along this mirror symmetry plane. Note that any incident polarization can be decomposed into two components that are parallel and perpendicular to this plane, respectively. Second, when a plane wave is normally incident on a nanostructure such that the field vector (or its components) are parallel to one mirror symmetry plane and perpendicular to another, the simulation domain can be cut to one quarter along these two perpendicular mirror symmetry planes [2]. In the case of a periodic unit cell, it is additionally required that such mirror planes bisecting the unit cell need to be parallel with the unit cell boundaries, i.e., the unit cell needs to be rectangular (not necessarily the primitive cell, however). The remaining boundaries opposite to the cut-planes are then assigned the same PMC or PEC boundary conditions. In the case of a single nanostructure, such a restriction does not apply and the remaining boundaries are typically terminated with perfectly matched layers (PMLs) to represent open boundaries extending to infinity.

An additional symmetry reduction applies in certain cases when considering the absorption of a normally incident plane wave. First, with a single nanostructure having a circular or a regular polygon cross-section, the absorption becomes independent of the incident field polarization. Second, similar polarization independence can be found for nanostructures arranged in a tetragonal or hexagonal lattice when the nanostructures have the corresponding symmetries or circular cross-section. In both cases, the polarization independence is simply due to there being two (or more) non-collinear, symmetry-equivalent directions to which any incident polarization can be decomposed. In general, however, single nanostructures and nanostructure arrays show varying scattering and absorption depending on the polarization and direction of the incident plane wave.

2.3. FEM Model Symmetry Reduction with Symmetry Modes

In the following, we will discuss how it is actually possible to similarly reduce the FEM model simulation domain to one half or one quarter with mirror symmetry planes of the geometry also in such cases when the electric field does not coincide with these geometry symmetries and hence does not fulfill the assigned PMC or PEC boundary conditions. This is achieved by expressing the electric field as a linear combination of symmetry modes which by themselves are solutions to Equation (1) and fulfill the boundary conditions of the reduced simulation domain (by virtue of the symmetry). Therefore, instead of solving the OP, we can solve SPs with each symmetry mode in the reduced simulation domain and combine the solutions from these to obtain the OP solution.

2.3.1. The Method

By selecting a suitable basis and using well-known results from group theory (the reduction formula and the projection operator method [10]), it is possible to find the symmetry modes and perform the decomposition. Such a procedure in the case of the ${\mathrm{C}}_{\mathrm{s}}$ and ${\mathrm{C}}_{2\mathrm{v}}$ point groups has already been carried out in detail elsewhere [17,18], so we will not include the derivations here. With the ${\mathrm{C}}_{2\mathrm{v}}$ point group, one can take as a basis the set $\{\widehat{\mathrm{E}}\mathbf{E}(x,y,z),{\widehat{\mathrm{C}}}_2\mathbf{E}(x,y,z),{\widehat{\sigma }}_{xz}\mathbf{E}(x,y,z),{\widehat{\sigma }}_{yz}\mathbf{E}(x,y,z)\}$, where the symmetry operations in ${\mathrm{C}}_{2\mathrm{v}}$ have been applied to the electric field solution $\mathbf{E}(x,y,z)$ in the geometry (the set of functions spanning a 4-dimensional function space). Clearly, $\mathbf{E}(x,y,z)$ can be trivially expressed on this basis, and performing any symmetry operation in ${\mathrm{C}}_{2\mathrm{v}}$ to any function in this basis yields one of the other functions in the same basis due to the closure of the group. With this starting point, the following un-normalized symmetry modes

$$\left\{\begin{array}{c}{\mathbf{E}}_1(x,y,z)=(\widehat{\mathrm{E}}+{\widehat{\mathrm{C}}}_2+{\widehat{\sigma }}_{xz}+{\widehat{\sigma }}_{yz})\mathbf{E}(x,y,z)\hfill \\ {\mathbf{E}}_2(x,y,z)=(\widehat{\mathrm{E}}+{\widehat{\mathrm{C}}}_2-{\widehat{\sigma }}_{xz}-{\widehat{\sigma }}_{yz})\mathbf{E}(x,y,z)\hfill \\ {\mathbf{E}}_3(x,y,z)=(\widehat{\mathrm{E}}-{\widehat{\mathrm{C}}}_2+{\widehat{\sigma }}_{xz}-{\widehat{\sigma }}_{yz})\mathbf{E}(x,y,z)\hfill \\ {\mathbf{E}}_4(x,y,z)=(\widehat{\mathrm{E}}-{\widehat{\mathrm{C}}}_2-{\widehat{\sigma }}_{xz}+{\widehat{\sigma }}_{yz})\mathbf{E}(x,y,z)\hfill \end{array}\right.$$

(2)

can be obtained (see ref. [18] for details) and the decomposition becomes

$$\mathbf{E}(x,y,z)=\frac{1}{4}({\mathbf{E}}_1(x,y,z)+{\mathbf{E}}_2(x,y,z)+{\mathbf{E}}_3(x,y,z)+{\mathbf{E}}_4(x,y,z)).$$

(3)

Similarly, with the ${\mathrm{C}}_{\mathrm{s}}$ point group, the following un-normalized symmetry modes

$$\left\{\begin{array}{c}{\mathbf{E}}_1(x,y,z)=(\widehat{\mathrm{E}}+{\widehat{\sigma }}_{xz})\mathbf{E}(x,y,z)\hfill \\ {\mathbf{E}}_2(x,y,z)=(\widehat{\mathrm{E}}-{\widehat{\sigma }}_{xz})\mathbf{E}(x,y,z)\hfill \end{array}\right.$$

(4)

can be obtained (see ref. [18] for details) and the decomposition becomes

$$\mathbf{E}(x,y,z)=\frac{1}{2}({\mathbf{E}}_1(x,y,z)+{\mathbf{E}}_2(x,y,z)).$$

(5)

For the FEM implementation, it then remains to determine which of the symmetry modes correspond to which combination of PMC and PEC boundary conditions for the mirror symmetry planes in the corresponding SPs. First, we consider the effect of the symmetry operations on $\mathbf{E}(x,y,z)$:

$$\left\{\begin{array}{c}\widehat{\mathrm{E}}\mathbf{E}(x,y,z)=({\mathbf{u}}_x{\mathbf{u}}_x+{\mathbf{u}}_y{\mathbf{u}}_y+{\mathbf{u}}_z{\mathbf{u}}_z)·\mathbf{E}(x,y,z)\hfill \\ {\widehat{\mathrm{C}}}_2\mathbf{E}(x,y,z)=(-{\mathbf{u}}_x{\mathbf{u}}_x-{\mathbf{u}}_y{\mathbf{u}}_y+{\mathbf{u}}_z{\mathbf{u}}_z)·\mathbf{E}(-x,-y,z)\hfill \\ {\widehat{\sigma }}_{xz}\mathbf{E}(x,y,z)=({\mathbf{u}}_x{\mathbf{u}}_x-{\mathbf{u}}_y{\mathbf{u}}_y+{\mathbf{u}}_z{\mathbf{u}}_z)·\mathbf{E}(x,-y,z)\hfill \\ {\widehat{\sigma }}_{yz}\mathbf{E}(x,y,z)=(-{\mathbf{u}}_x{\mathbf{u}}_x+{\mathbf{u}}_y{\mathbf{u}}_y+{\mathbf{u}}_z{\mathbf{u}}_z)·\mathbf{E}(-x,y,z)\hfill \end{array}\right.,$$

(6)

where ${\mathbf{u}}_x{\mathbf{u}}_x$, ${\mathbf{u}}_y{\mathbf{u}}_y$, and ${\mathbf{u}}_z{\mathbf{u}}_z$ are the unit vector dyads in the Cartesian coordinate system (this notation is adopted here to keep the expressions concise). Applying Equation (6) to the ${\mathrm{C}}_{2\mathrm{v}}$ symmetry modes of Equation (2) allows us to find the relation between these modes and the boundary conditions on the $xz$ and $yz$ planes as summarized in

Table 1. Finding the $xz$ and $yz$ plane boundary conditions for each symmetry mode corresponding to the ${\mathrm{C}}_{2\mathrm{v}}$ point group.

Mode	$E(\mathit{x},0,\mathit{z})\times u_{\mathit{y}}$	$E(\mathit{x},0,\mathit{z})·u_{\mathit{y}}$	$E(0,\mathit{y},\mathit{z})\times u_{\mathit{x}}$	$E(0,\mathit{y},\mathit{z})·u_{\mathit{x}}$	$\mathit{xz}$ Plane	$\mathit{yz}$ Plane
${\mathbf{E}}_1(x,y,z)$	$\ne 0$	$=0$	$\ne 0$	$=0$	PMC	PMC
${\mathbf{E}}_2(x,y,z)$	$=0$	$\ne 0$	$=0$	$\ne 0$	PEC	PEC
${\mathbf{E}}_3(x,y,z)$	$\ne 0$	$=0$	$=0$	$\ne 0$	PMC	PEC
${\mathbf{E}}_4(x,y,z)$	$=0$	$\ne 0$	$\ne 0$	$=0$	PEC	PMC

. Similarly, applying Equation (6) to the ${\mathrm{C}}_{\mathrm{s}}$ symmetry modes of Equation (4) allows us to find the relation between these modes and the boundary conditions on the $xz$ plane as summarized in

Table 2. Finding the $xz$ plane boundary condition for each symmetry mode corresponding to the ${\mathrm{C}}_{\mathrm{s}}$ point group.

Mode	$\mathbf{E}(\mathit{x},0,\mathit{z})\times {\mathbf{u}}_{\mathit{y}}$	$\mathbf{E}(\mathit{x},0,\mathit{z})·{\mathbf{u}}_{\mathit{y}}$	$\mathit{xz}$ Plane
${\mathbf{E}}_1(x,y,z)$	$\ne 0$	$=0$	PMC
${\mathbf{E}}_2(x,y,z)$	$=0$	$\ne 0$	PEC

.

We use the background-field-scattered-field formulation in our FEM models, where the OP is defined by the geometry, refractive indices, boundary conditions, and a background field solution in the absence of the nanostructures. The electric field solution can therefore be expressed as $\mathbf{E}(x,y,z)={\mathbf{E}}_{\mathrm{bg}}(x,y,z)+{\mathbf{E}}_{\mathrm{s}}(x,y,z)$, where ${\mathbf{E}}_{\mathrm{s}}(x,y,z)$ is the scattered field solved in the FEM process. This means that, due to the linearity of Maxwell’s equations, we can actually apply the symmetry mode decomposition (Equation (3) or Equation (5)) to the known background field solution ${\mathbf{E}}_{\mathrm{bg}}(x,y,z)$ when we set up the SPs. Therefore, as long as we know the OP background field solution, we can apply the symmetry reduction and obtain fully defined SPs. Note that we have assumed here that the symmetry modes obtained for the background field as described above are solutions to Equation (1) themselves. This assumption is, indeed, true when the incident field is a plane wave or, more generally, can be expressed as a plane wave series expansion, since applying the considered symmetry operators (Equation (6)) to plane wave components satisfying Equation (1) produces other plane waves that still satisfy Equation (1).

When the SPs have been solved, the OP solution in the reduced simulation domain is obtained with either Equation (3) or Equation (5) in the case of ${\mathrm{C}}_{2\mathrm{v}}$ or ${\mathrm{C}}_{\mathrm{s}}$ symmetry, respectively. However, the solution is needed in the full OP simulation domain. Since the solutions of each SP show different combinations of symmetric and antisymmetric fields with respect to the mirror planes (symmetric with PMC and antisymmetric with PEC boundary condition, see Table 1 and Table 2), care must be taken with the signs when summing them to construct the OP solution in the full domain. If we take the reduced domain to be the first quadrant ($x\ge 0,y\ge 0$) with the ${\mathrm{C}}_{2\mathrm{v}}$ symmetry or half ($y\ge 0$) with the ${\mathrm{C}}_{\mathrm{s}}$ symmetry, the full domain OP solution can be expressed as

$$\mathbf{E}(x,y,z)=\left\{\begin{array}{cc}\frac{1}{4}\left({\mathbf{E}}_1(x,y,z)+{\mathbf{E}}_2(x,y,z)+{\mathbf{E}}_3(x,y,z)+{\mathbf{E}}_4(x,y,z)\right)\hfill & ,x\ge 0,y\ge 0\hfill \\ \frac{1}{4}\left({\widehat{\sigma }}_{yz}{\mathbf{E}}_1(x,y,z)-{\widehat{\sigma }}_{yz}{\mathbf{E}}_2(x,y,z)-{\widehat{\sigma }}_{yz}{\mathbf{E}}_3(x,y,z)+{\widehat{\sigma }}_{yz}{\mathbf{E}}_4(x,y,z)\right)\hfill & ,x<0,y\ge 0\hfill \\ \frac{1}{4}\left({\widehat{\mathrm{C}}}_2{\mathbf{E}}_1(x,y,z)+{\widehat{\mathrm{C}}}_2{\mathbf{E}}_2(x,y,z)-{\widehat{\mathrm{C}}}_2{\mathbf{E}}_3(x,y,z)-{\widehat{\mathrm{C}}}_2{\mathbf{E}}_4(x,y,z)\right)\hfill & ,x<0,y<0\hfill \\ \frac{1}{4}\left({\widehat{\sigma }}_{xz}{\mathbf{E}}_1(x,y,z)-{\widehat{\sigma }}_{xz}{\mathbf{E}}_2(x,y,z)+{\widehat{\sigma }}_{xz}{\mathbf{E}}_3(x,y,z)-{\widehat{\sigma }}_{xz}{\mathbf{E}}_4(x,y,z)\right)\hfill & ,x\ge 0,y<0\hfill \end{array}\right.$$

(7)

or as

$$\mathbf{E}(x,y,z)=\left\{\begin{array}{cc}\frac{1}{2}\left({\mathbf{E}}_1(x,y,z)+{\mathbf{E}}_2(x,y,z)\right)\hfill & ,y\ge 0\hfill \\ \frac{1}{2}\left({\widehat{\sigma }}_{xz}{\mathbf{E}}_1(x,y,z)-{\widehat{\sigma }}_{xz}{\mathbf{E}}_2(x,y,z)\right)\hfill & ,y<0\hfill \end{array}\right.,$$

(8)

respectively. Note that these expressions (and subsequent computations) for each quartile or half can actually still be evaluated in the reduced domain mesh, i.e., it is not necessary to construct the full domain mesh to process the OP solution (depending on the specifics of the FEM implementation, this may or may not simplify the postprocessing).

2.3.2. Plane Wave Incidence

We restrict our discussion here to models with an incident plane wave field. As illustrated in Figure 3, we define the plane wave incidence with three angles: polar angle ${\theta }_{\mathrm{i}}$ (from z-axis), azimuth angle ${\varphi }_{\mathrm{i}}$ (from x-axis), and electric field polarization rotation angle ${\psi }_{\mathrm{i}}$ (with ${\psi }_{\mathrm{i}}=0$ and ${\psi }_{\mathrm{i}}=\pi /2$ corresponding to p- and s-polarization, i.e., parallel and perpendicular to the plane of incidence, respectively). We also define the normal incidence as ${\theta }_{\mathrm{i}}=0$ and ${\varphi }_{\mathrm{i}}=0$ such that the polarization rotation is always defined solely by ${\psi }_{\mathrm{i}}$. The incident electric field can then be expressed as

$${\mathbf{E}}_{\mathrm{i}}(x,y,z)={\mathbf{E}}_{\mathrm{i}}{e}^{-i(k_{\mathrm{i}x}x+k_{\mathrm{i}y}y+k_{\mathrm{i}z}z)},$$

(9)

with the wave vector

$${\mathbf{k}}_{\mathrm{i}}=-n_1{k}_0\left[\begin{array}{c}sin\left({\theta }_{\mathrm{i}}\right)cos\left({\varphi }_{\mathrm{i}}\right)\\ sin\left({\theta }_{\mathrm{i}}\right)sin\left({\varphi }_{\mathrm{i}}\right)\\ cos\left({\theta }_{\mathrm{i}}\right)\end{array}\right]$$

(10)

and electric field vector

$${\mathbf{E}}_{\mathrm{i}}=E_0\left(cos\left({\psi }_{\mathrm{i}}\right)\left[\begin{array}{c}cos\left({\theta }_{\mathrm{i}}\right)cos\left({\varphi }_{\mathrm{i}}\right)\\ cos\left({\theta }_{\mathrm{i}}\right)sin\left({\varphi }_{\mathrm{i}}\right)\\ -sin\left({\theta }_{\mathrm{i}}\right)\end{array}\right]+sin\left({\psi }_{\mathrm{i}}\right)\left[\begin{array}{c}-sin\left({\varphi }_{\mathrm{i}}\right)\\ cos\left({\varphi }_{\mathrm{i}}\right)\\ 0\end{array}\right]\right),$$

(11)

where $n_1$ is the refractive index of the incidence medium, $k_0$ is the wave number in free space, $E_0$ is the electric field amplitude, and the column vectors represent the x-, y-, and z-components.

In a general stratified background medium, the incident plane wave results in a field solution of the form

$${\mathbf{E}}_{\mathrm{bg}}=\left[\begin{array}{c}{E}_x\left(z\right)\\ E_y\left(z\right)\\ E_z\left(z\right)\end{array}\right]e^{-i(k_{\mathrm{i}x}x+k_{\mathrm{i}y}y)},$$

(12)

where the z-dependence can be solved with the transfer matrix method (see e.g., ref. [19] and the references therein). As with plane wave components in a homogeneous medium, applying the considered symmetry operators (Equation (6)) to the stratified medium background field solution (Equation (12)) produce fields that still satisfy Equation (1), and the obtained background field symmetry modes are hence solutions to Equation (1) also in this case. In our numerical examples, we consider the simple case of a single interface between the superstrate and substrate (with refractive indices $n_1$ and $n_2$, respectively) resulting in a single reflected and transmitted plane wave component with the well-known Fresnel reflection and transmission coefficients.

2.3.3. Periodic Nanostructure Array

In the case of a single nanostructure, the above considerations apply with a general plane wave incidence, but additional restrictions arise when we try to apply this method to symmetry reduce a unit cell of a periodic nanostructure array. With the unit cell, we need to consider not only the $xz$ and $yz$ mirror plane symmetries but also the discrete translational symmetry of the lattice. We take the unit cell to be centered at the origin and the lattice period to be $p_x$ in x-direction and $p_y$ in y-direction. In the case of ${\mathrm{C}}_{2\mathrm{v}}$ symmetry, in order to use the previous symmetry reduction of the simulation domain to the first quadrant $\{0\le x\le p_x/2,0\le y\le p_y/2\}$ and apply the same boundary conditions on opposite boundaries, we find that the plane wave incidence must fulfill the conditions

$$\left\{\begin{array}{cc}{k}_{\mathrm{i}x}=-\frac{2\pi }{p_x}{m}_x\hfill & ,m_x\in \mathbb{Z}\hfill \\ k_{\mathrm{i}y}=-\frac{2\pi }{p_y}{m}_y\hfill & ,m_y\in \mathbb{Z}\hfill \end{array}\right.$$

(13)

for the boundary conditions on the $x=p_x/2$ and $y=p_y/2$ boundaries also to be satisfied. These incidence conditions result since, for a plane wave background, the symmetry modes consist of plane wave components which then need to have zero or multiple of $2\pi$ phase difference on these boundaries in order for the resultant sum field to satisfy the PMC or PEC boundary conditions. In the case of ${\mathrm{C}}_{\mathrm{s}}$ symmetry with the simulation domain reduced to the half $\{-p_x/2\le x\le p_x/2,0\le y\le p_y/2\}$, Floquet boundary conditions can still be applied on the $x=\pm p_x/2$ boundaries and only the condition $k_{\mathrm{i}y}=-\frac{2\pi }{p_y}{m}_y$ is required for the PMC or PEC boundary condition on the $y=0$ and $y=p_y/2$ boundaries to be satisfied. It is also important to note here that, for lattices where the primitive cell is not rectangular, the smallest rectangular unit cell required for this symmetry reduction method is larger than the primitive cell and the obtained benefit is reduced, especially since using the primitive cell and Floquet boundary conditions allows for arbitrary incidence. For example, in a hexagonal lattice (see Figure 2c), the smallest rectangular unit cell is twice the size of the primitive cell and the symmetry reduction to one quarter reduces the simulation domain size by only a factor of 2 compared to the primitive cell.

We note that the incidence conditions of Equation (13) actually appear in the context of diffraction gratings. Following from the periodicity of the geometry and the incident plane wave, the wave vector components of the diffracted plane waves (in either the superstrate or the substrate) of diffraction order $({\nu }_x,{\nu }_y)$ (with ${\nu }_x$ and ${\nu }_y$ integers) need to satisfy the conditions

$$\left\{\begin{array}{c}{k}_x^{\left({\nu }_x\right)}=k_{\mathrm{i}x}+\frac{2\pi }{p_x}{\nu }_x\hfill \\ k_y^{\left({\nu }_y\right)}=k_{\mathrm{i}y}+\frac{2\pi }{p_y}{\nu }_y\hfill \end{array}\right.$$

(14)

and the incidence condition of Equation (13) can therefore be interpreted as the so-called Littrow configuration of order $(2m_x,2m_y)$ (i.e., a configuration where the reflected diffraction order propagates back to the incidence direction: ${\mathbf{k}}^{({\nu }_x,{\nu }_y)}=-{\mathbf{k}}_{\mathrm{i}}$). Furthermore, finding such a restriction is not unique to our FEM approach since similar Littrow configuration incidence condition requirements were found in ref. [18] with ${\mathrm{C}}_{2\mathrm{v}}$ symmetry reduction for FMM. Indeed, this requirement is rather a result of the symmetry itself than of the chosen simulation method.

2.3.4. Special Incidence Cases

The symmetry reduction method outlined above has some special plane wave incidence cases including the previously discussed well-known ones of normal incidence for ${\mathrm{C}}_{2\mathrm{v}}$ symmetry and the plane of incidence coinciding with the mirror symmetry plane for ${\mathrm{C}}_{\mathrm{s}}$ symmetry. With normal incidence and ${\mathrm{C}}_{2\mathrm{v}}$ symmetry, ${\mathbf{E}}_1(x,y,z)={\mathbf{E}}_2(x,y,z)=0$, ${\mathbf{E}}_3(x,y,z)=4{\mathbf{u}}_x{\mathbf{u}}_x·\mathbf{E}(x,y,z)$, and ${\mathbf{E}}_4(x,y,z)=4{\mathbf{u}}_y{\mathbf{u}}_y·\mathbf{E}(x,y,z)$, reducing the background field decomposition and SPs to just consider the two perpendicular polarization components separately. With incidence along the symmetry plane (${\varphi }_{\mathrm{i}}=0$) and ${\mathrm{C}}_{\mathrm{s}}$ symmetry, ${\mathbf{E}}_1(x,y,z)$ and ${\mathbf{E}}_2(x,y,z)$ become the p- and s-polarized component, respectively, reducing the background field decomposition and SPs to just consider the two polarization components separately.

However, in a third special case of symmetry plane incidence (${\varphi }_{\mathrm{i}}=0$ or ${\varphi }_{\mathrm{i}}=\pi /2$) and ${\mathrm{C}}_{2\mathrm{v}}$ symmetry, in general, all four symmetry modes are needed but, with the incident field either fully p- or s-polarized (${\psi }_{\mathrm{i}}=0$ or ${\psi }_{\mathrm{i}}=\pi /2$), only two of the four symmetry modes are nonzero, reducing the number of SPs to consider to two as listed in

Table 3. SP symmetry modes in the special plane wave incidence cases where the plane of incidence coincides with either the $xz$ or $yz$ mirror symmetry plane (with ${\mathrm{C}}_{2\mathrm{v}}$ symmetry).

Incidence Plane	Polarization	Nonzero Symmetry Modes
$xz$	p	${\mathbf{E}}_1$ and ${\mathbf{E}}_3$
$xz$	s	${\mathbf{E}}_2$ and ${\mathbf{E}}_4$
$yz$	p	${\mathbf{E}}_1$ and ${\mathbf{E}}_4$
$yz$	s	${\mathbf{E}}_2$ and ${\mathbf{E}}_3$

. It is rather straightforward to set up the symmetry reduced models to identify these special incidence cases and only solve the needed SPs.

2.3.5. Expected Performance Increase

From computational cost point of view, we on the one hand reduce the simulation domain by a factor of four or two with ${\mathrm{C}}_{2\mathrm{v}}$ or ${\mathrm{C}}_{\mathrm{s}}$ symmetry, respectively, but on the other hand we increase the number of problems to solve from the one OP to four or two SPs, respectively. The degrees of freedom (DOF) and hence the FEM system matrix size are expected to scale approximately linearly with the simulation domain size. The full FEM system matrix size actually scales as the square of DOF, but this matrix is sparse since the Equation (1) is local, connecting only neighboring mesh elements to each other. Therefore, the number of non-zero matrix elements is proportional to DOF and so is the random access memory (RAM) use with sparse solvers that store only the non-zero elements. In our previous numerical simulation benchmarking study [20] we found the maximum RAM usage to scale approximately linearly with DOF and the solver run time to scale superlinearly. However, for small DOF, overhead starts to be significant and neither the maximum RAM usage nor run time will scale linearly. Additionally, the solver run time, in general, depends heavily on the number and type of CPU cores utilized.

It is also worth taking into consideration how the simulation domain boundaries affect the DOF and RAM use. PMLs in non-periodic models are not boundary conditions but additional domains (within which the fields are absorbed without reflection) and hence increase the simulation domain size and DOF. Floquet periodic boundary conditions in periodic models instead connect the field components on opposite boundaries making the system matrix more dense without increasing DOF. Finally, the PMC and PEC boundary conditions simply set some field components to zero on the boundary, effectively decreasing the DOF and hence result in the lowest RAM use compared to using other boundary conditions.

Overall, it is to be expected that our symmetry reduction method, with non-periodic models, should reduce the memory requirement approximately by a factor of 4 or 2 for ${\mathrm{C}}_{2\mathrm{v}}$ or ${\mathrm{C}}_{\mathrm{s}}$ symmetry, respectively, and to possibly also yield reduced total solver run time, especially with high DOF. With periodic models, our symmetry reduction method should additionally reduce the RAM use and total solver time due to changing Floquet boundary conditions to PMC or PEC. Furthermore, in the special incidence cases mentioned above, the number of SPs to solve is reduced leading to reduced total solver run time.

3. Numerical Examples

We demonstrate our FEM model symmetry reduction with symmetry modes using two numerical examples, implemented with COMSOL Multiphysics^® 5.6 software with the Wave Optics Module. First, we consider models of a single hexagonal nanowire embedded in a semi-infinite superstrate on top of a a semi-infinite substrate for obtaining the nanowire absorption cross-section as a function of the plane wave incidence direction and polarization (as opposed to the common practice of modeling at just the normal incidence for convenience). Such models could be used e.g., in the analysis or design of a single-nanowire solar cell [3]. Second, we consider models also with semi-infinite superstrate and substrate but with a periodic array of hexagonal nanowires arranged in an orthorhombic (rectangular) lattice. In this case, the primitive cell is rectangular, and we therefore get the full benefit of using the symmetry reduction (as opposed to a non-rectangular primitive cell).

Due to the restriction on the plane wave incidence angle (Equation (13)), this second example with symmetry reduced unit cell for the array finds somewhat limited applicability. One application could be to use such an array as a grating and use similar models to compute the diffraction order efficiencies as a function of the nanowire dimensions (diameter and length) when the incidence is kept fixed at one of the $(2m_x,2m_y)$th Littrow configurations (for chosen wavelength and period). Alternatively, the array could act as a photodetector [5] and the models be used to design the nanowire dimensions to maximize absorption at the desired wavelengths (again, with a fixed specific incidence that respects the restrictions). We lean towards the latter application in our example but also demonstrate an additional aspect of the symmetry reduction.

Therefore, in our periodic nanowire array example, we consider the incidence of two coherent antiphase plane waves at mirror incidence angles (${\theta }_{\mathrm{i}2}={\theta }_{\mathrm{i}1}$, ${\varphi }_{\mathrm{i}1}=0$, ${\varphi }_{\mathrm{i}2}=\pi$, and ${\psi }_{\mathrm{i}2}=\pi -{\psi }_{\mathrm{i}1}$) leading to polarization selective absorption in the nanowires (centered at the origin of the unit cell) with constructive interference for the p-polarization x-components and destructive interference for the s-polarization y-components (also destructive interference for the p-polarization z-components). We note that, with this incidence configuration and either p- or s-polarization, the background field is actually directly one of the symmetry modes: ${\mathbf{E}}_3(x,y,z)$ for p- and ${\mathbf{E}}_2(x,y,z)$ for s-polarization, respectively. Therefore, this is a simple example of such a case where the symmetry modes themselves could be of interest. We further note that exploiting such coherent absorption of multiple incident plane waves has wider interest [21] and similar opportunities for symmetry reduction of models could arise in this context. With these models, we explore the array absorption as a function of the periods $p_x$ and $p_y$, while maintaining the $(2,0)$ and $(-2,0)$ Littrow configuration for the two incident plane waves, i.e., $k_{\mathrm{i}x}=\pm 2\pi /p_x$ and $k_{\mathrm{i}y}=0$, meaning that $p_y$ does not put restriction on the incidence angle.

3.1. Models

In the first single nanowire absorption cross-section model, we select the substrate material as Si, the nanowire material as GaAs, the nanowire diameter as $D=380$ nm, and the nanowire length as $L=2500$ nm (viz. ref. [3]). The wavelength (in free space) of the incident plane wave is chosen as ${\lambda }_0=800$ nm, at which we use the refractive index $n=3.6750-0.0054113i$ for Si [22] and $n=3.6035-0.091786i$ for GaAs [23]. For simplicity, we consider the superstrate medium to be a polymer with a refractive index of $n=1.4$ instead of a more realistic stack of air, transparent conductive oxide, and polymer (viz. ref [3]). As already mentioned, with this simplification, the general background field expression of Equation (12) reduces to a total of three plane wave components (incident and reflected plane wave in the superstrate and transmitted plane wave in the substrate). We select the simulation domain width (square cross-section), superstrate height, and substrate height such that they are sufficient to not result in numerical artifacts in the field solution (approximately one wavelength space between the nanowire and the simulation domain boundaries). We then cut the simulation domain with the $xz$ and $yz$ planes to one quarter with $x\ge 0$ and $y\ge 0$. All the remaining original boundaries are terminated with PMLs, while for each SP, the background field is set as one of the symmetry modes of Equation (2) and the cut plane boundaries are assigned either PMC or PEC boundary condition according to Table 1. The simulation domain is meshed with a free tetrahedral mesh using the Fine preset settings in the software and adjusting the maximum mesh element size to check convergence of the results with mesh refinement. The PMLs are meshed with a swept mesh of six elements across. We selected the tightest maximum mesh element size as ${\lambda }_0/\left(12n\right)$, where ${\lambda }_0$ is the free space wavelength and n is the refractive index (real part) in each domain, and found sufficient convergence in the obtained absorption cross-section values with maximum mesh element size of ${\lambda }_0/\left(3n\right)$. The electric field in the model is solved using a Wavelength Domain study at ${\lambda }_0=800$ nm and the direct MUMPS solver with the default settings, except for turning off the out-of-core functionality to avoid swapping data from RAM to hard drive storage if the system matrix becomes too large, which would lead to erroneous extracted values of maximum RAM use. We opted to use a direct solver instead of an iterative one since, although iterative solvers tend to use less RAM, direct solvers are more robust and reliable. For comparison, we also consider a second model without the symmetry reduction having the OP background field and all boundaries terminated with PMLs. In the following, we refer to these two models as “Sym. red. 1/4” and “No sym. red.”, respectively.

In the first periodic nanowire array absorption model, we select the substrate as Si, the nanowire material as In${}_{0.52}$Ga${}_{0.48}$As, the nanowire diameter as $D=260$ nm, and the nanowire length as $L=1400$ nm (viz. ref. [5], but we choose a thicker nanowire for increased absorption, although not optimized [8]). The wavelength (in free space) of the incident plane wave is chosen as the popular telecom wavelength ${\lambda }_0=1550$ nm, at which we use the refractive index $n=3.4757$ for Si [24] and $n=3.5307-0.076399i$ for In${}_{0.52}$Ga${}_{0.48}$As [23]. Again, for simplicity, we consider the superstrate medium to be a polymer with a refractive index of $n=1.4$ instead of a more realistic stack of air, transparent conductive oxide, and polymer (viz. ref [5]). We select the superstrate and substrate heights such that they are sufficient to not result in numerical artifacts in the field solution (approximately one wavelength space between the nanowire and the simulation domain boundaries above and below). We then cut the simulation domain with the $xz$ and $yz$ planes to one quarter with $x\ge 0$ and $y\ge 0$. The top and bottom boundaries are terminated with PMLs, and the background field with the two incident plane waves is, as mentioned, the symmetry mode ${\mathbf{E}}_3(x,y,z)$ for p- and ${\mathbf{E}}_2(x,y,z)$ for s-polarization (Equation 2 and the cut plane boundaries are assigned either PMC or PEC boundary condition according to Table 1). The simulation domain is meshed with a free tetrahedral mesh using the Fine preset settings and adjusting the maximum mesh element size to check convergence of the results with mesh refinement. The PMLs are meshed with a swept mesh of six elements across. We selected the tightest maximum mesh element size as ${\lambda }_0/\left(12n\right)$, where ${\lambda }_0$ is the free space wavelength and n is the refractive index (real part) in each domain, and found sufficient convergence in the obtained absorption values with maximum mesh element size of ${\lambda }_0/\left(6n\right)$, except for small $p_y$. The electric fields in the model are solved using a Wavelength Domain study at ${\lambda }_0=1550$ nm and the direct MUMPS solver with the default settings, except for turning off the out-of-core functionality. For the periodic models, iterative solvers are not even a practical option since they, in general, tend to not work well with Floquet boundary conditions. Therefore, we use the selected direct MUMPS solver throughout. In the following, this first periodic model is referred to as “APMI sym. red. 1/4”.

For comparison, we also consider four additional models: one with the same background field but without the symmetry reduction and three single plane wave background field ($k_{\mathrm{i}x}=-2\pi /p_x$ and $k_{\mathrm{i}y}=0$) models: without symmetry reduction, with symmetry reduction to one half, and with symmetry reduction to one quarter, respectively. In the following, we refer to these four models as “APMI no sym. red.”, “SI no sym. red.”, “SI sym. red. 1/2”, and “SI sym. red. 1/4”, respectively. The model “APMI no sym. red.” has the original simulation domain with Floquet boundaries at the sides. Note that, since the wave vector x-components of the two incident plane waves differ only in sign and the y-components are zero, ${\mathbf{k}}_{\mathrm{F}}=k_{\mathrm{i}x}{\mathbf{u}}_x$ and the Floquet condition on the boundaries $y=\pm p_y/2$ reduces to field continuity. The only difference between the models “SI no sym. red.” and “APMI no sym. red.” is the background field. These two models without symmetry reduction serve as the reference point for the symmetry reduced models. The “SI sym. red. 1/2” model uses the symmetry modes of Equation (4) which simply reduce to the p- and s-polarized component of the background field. Boundary conditions for the $y=0$ and $y=p_y/2$ planes are determined from Table 2 while the $x=-p_x/2$ and $x=p_x/2$ planes have the Floquet boundary condition. This model actually corresponds to the aforementioned conventional mirror plane incidence symmetry reduction case where the incident light is along the mirror plane. Finally, background field in the “SI sym. red. 1/4” corresponds to the aforementioned special case of the symmetry modes of Equation (2), where ${\mathbf{E}}_2(x,y,z)={\mathbf{E}}_4(x,y,z)=0$ for p-polarization and ${\mathbf{E}}_1(x,y,z)={\mathbf{E}}_3(x,y,z)=0$ for s-polarization, reducing the SPs to two for each polarization. Otherwise, the four models are set up the same as the “APMI sym. red. 1/4” model. For clarity, we have summarized the SPs solved with each model in

Table 4. Summary of the SPs to solve in each periodic nanowire array absorption model depending on the polarization of the OP incident plane wave (see the main text for the definitions of the models and the symmetry modes ${\mathbf{E}}_j$ for ${\mathrm{C}}_{\mathrm{s}}$ and ${\mathrm{C}}_{2\mathrm{v}}$ symmetry).

Model	P-Polarized OP	S-Polarized OP	General OP
APMI sym. red. 1/4	${\mathbf{E}}_3$ (${\mathrm{C}}_{2\mathrm{v}}$)	${\mathbf{E}}_2$ (${\mathrm{C}}_{2\mathrm{v}}$)	${\mathbf{E}}_3$ and ${\mathbf{E}}_2$ (${\mathrm{C}}_{2\mathrm{v}}$)
APMI no sym. red.	OP	OP	OP
SI no sym. red.	OP	OP	OP
SI sym. red. 1/2	${\mathbf{E}}_1$ (${\mathrm{C}}_{\mathrm{s}}$)	${\mathbf{E}}_2$ (${\mathrm{C}}_{\mathrm{s}}$)	${\mathbf{E}}_1$ and ${\mathbf{E}}_2$ (${\mathrm{C}}_{\mathrm{s}}$)
SI sym. red. 1/4	${\mathbf{E}}_1$ and ${\mathbf{E}}_3$ (${\mathrm{C}}_{2\mathrm{v}}$)	${\mathbf{E}}_2$ and ${\mathbf{E}}_4$ (${\mathrm{C}}_{2\mathrm{v}}$)	${\mathbf{E}}_1$, ${\mathbf{E}}_2$, ${\mathbf{E}}_3$, and ${\mathbf{E}}_4$ (${\mathrm{C}}_{2\mathrm{v}}$)

.

With the field solutions from the solved models, we compute the single nanowire absorption cross-section and nanowire array absorption results. The nanowire absorption cross-section is computed as $A_{\mathrm{cs}}=P_{\mathrm{nw}}/I_{\mathrm{i}}$, where $P_{\mathrm{nw}}$ is the power absorbed in the nanowire ($\Re \left\{\mathbf{J}\right(x,y,z)·\mathbf{E}(x,y,z\left)\right\}$ integrated over the nanowire volume, where $\mathbf{J}(x,y,z)$ is the induced current density) and $I_{\mathrm{i}}$ is the incident intensity ($I_{\mathrm{i}}=n_1/\left(2{\eta }_0\right)E_0^2$, where ${\eta }_0$ is the free space impedance). The nanowire array absorption is computed as $A=P_{\mathrm{nw}}/P_{\mathrm{i}}$, where $P_{\mathrm{i}}$ is the incident power into the unit cell and $P_{\mathrm{nw}}$ is the power absorbed in the nanowire in the unit cell. Note that there is no power flow between the unit cells, so the incident power in one unit cell is given by integrating the incident field Poynting vector z-component over the cross-section area ($P_{\mathrm{i}}=p_x{p}_y{I}_{\mathrm{i}}cos\left({\theta }_{\mathrm{i}}\right)$). In the “APMI” models, the incident power is doubled with two plane waves incident instead of one (the interference does not alter the total power, only its distribution). We perform the computation of $P_{\mathrm{nw}}$ in the symmetry reduced domain with Equation (7) or Equation (8) and sum the results to obtain the value in the full domain.

3.2. Simulation Results

With the two single nanowire absorption cross-section models “Sym. red. 1/4” and “No sym. red.”, we considered both p- and s-polarization and swept the polar incidence angle ${\theta }_{\mathrm{i}}$ from $0^{\circ }$ to ${88}^{\circ }$ with $2^{\circ }$ step and the azimuth incidence angle ${\varphi }_{\mathrm{i}}$ from $0^{\circ }$ to ${30}^{\circ }$ with $1^{\circ }$ step (due to our definition of normal incidence, ${\varphi }_{\mathrm{i}}=0^{\circ }$ when ${\theta }_{\mathrm{i}}=0^{\circ }$ in the sweep). These ranges for the incidence angles covered all unique incidences from the upper half-space, due to the hexagonal nanowire cross-section (${\mathrm{C}}_{6\mathrm{v}}$ symmetry). Figure 4 shows the nanowire absorption cross-section for ${\varphi }_{\mathrm{i}}={30}^{\circ }$ obtained from the simulations and the maximum DOF, total solver time and maximum RAM used for the whole sweep. These reported results were obtained with a local desktop computer Dell Precision 3640 MT with Intel Core i7-10700 8-core CPU and 128 GB RAM. As seen in Figure 4a, the good agreement in the absorption cross-section values between the two models indicates that the used symmetry reduction was set up correctly.

The results shown in Figure 4b were in line with our expectations. Since the symmetry reduced domain in the “Sym. red. 1/4” model was a factor of 4 smaller than the full domain in the “No sym. red.” model, the DOF was also approximately a factor of 4 smaller and so was the maximum RAM use, accordingly. Interestingly, a superlinear dependence of the solver time on the system matrix size showed up in the results yielding also a much smaller total solver time for the “Sym. red. 1/4” model. We did take into account the special incidence angle cases in the sweep where only part of the SPs needed to be solved, but this concerned only the normal incidence and ${\varphi }_{\mathrm{i}}=0^{\circ }$ (which constituted only approximately 1/30 of the steps in the sweep) and was hence not the reason for the observed reduced total solver time with symmetry reduction. We noticed that the overhead from setting up the solver was increased for the “Sym. red. 1/4” model when several SPs were solved instead of the one OP. However, as long as the overhead time was small compared to the actual solver time, as is likely to be the case with tight mesh for good convergence, this was not an issue for obtaining performance increase with the symmetry reduction.

With the periodic nanowire array absorption models “APMI sym. red. 1/4”, “APMI no sym. red.”, “SI no sym. red.”, “SI sym. red. 1/2”, and “SI sym. red. 1/4”, we considered both p- and s-polarization and swept the period $p_x$ such that, with the (2,0) Littrow configuration, the corresponding polar incidence angle ${\theta }_{\mathrm{i}}$ was swept from ${10}^{\circ }$ to ${60}^{\circ }$ with $0.5^{\circ }$ step. We first also considered different $p_y$ periods from 300 nm (less than one diameter D space between the nanowires) to 1107 nm, corresponding to ${\lambda }_0/n_1$, with coarser steps and meshing and then selected $p_y=950$ nm for this example due to good convergence in the absorption results with the maximum mesh element size of ${\lambda }_0/\left(6n\right)$. Figure 5 shows the periodic nanowire array absorption obtained from the simulations with the models “APMI sym. red. 1/4” and “APMI no sym. red.” and the maximum DOF, total solver time, and maximum RAM used for the sweep. Figure 6 shows the periodic nanowire array absorption obtained from the simulations with the models “SI no sym. red.”, “SI sym. red. 1/2”, and “SI sym. red. 1/4” and the maximum DOF, total solver time, and maximum RAM used for the $p_x$ (or equivalently ${\theta }_{\mathrm{i}}$) sweep. These reported results were also obtained with the local desktop computer Dell Precision 3640 MT with Intel Core i7-10700 8-core CPU and 128 GB RAM. The good agreement in the absorption values between the models in both sets (as seen in Figure 5a and Figure 6a) indicated that the used symmetry reductions were set up correctly. Furthermore, the expected difference in the polarization dependence of the absorption between the case of two antiphase plane waves at mirror incidence and the case of a single plane wave incidence was also demonstrated.

Interestingly, the results shown in Figure 5b and Figure 6b were much better than what could be expected based on the DOF reduction alone: although the DOF scaled as expected, the maximum RAM use and total solver time with all the symmetry reduced models were much smaller than that would entail. With the “APMI sym. red. 1/4” model the maximum RAM use was approximately a factor of six smaller and the total solver time was approximately a factor of 10 smaller than with the “APMI no sym. red.” model. With the “SI sym. red. 1/2” model the maximum RAM use was approximately a factor of three smaller and the total solver time was approximately a factor of five smaller than with the “SI no sym. red.” model. Finally, with the “SI sym. red. 1/4” model the maximum RAM use was approximately a factor of nine smaller and the total solver time was approximately a factor of seven smaller than with the “SI no sym. red.” model. We attribute this result to the aforementioned effect of Floquet (or field continuity) boundary conditions making the system matrix more dense than the simpler PMC or PEC boundary conditions. Indeed, the symmetry reduced periodic array models having the additional benefit of not only making the system matrix smaller but also less dense, clearly could exhibit strongly reduced computational cost.

4. Discussion

In summary, we have shown how to symmetry reduce an FEM optics model with nanostructures on top of a substrate and incident plane waves as the excitation by using the $xz$ and $yz$ mirror symmetry planes of the geometry and the associated symmetry modes to decompose the OP to SPs with reduced simulation domain size and computational cost. This symmetry reduction method can be applied when the model geometry either belongs to the ${\mathrm{C}}_{\mathrm{s}}$ or ${\mathrm{C}}_{2\mathrm{v}}$ point group or when it belongs to a point group that has one or both of these groups as a subgroup. It is then possible to solve the OP in the background-field-scattered-field formulation by taking the OP background field that does not need to coincide with the symmetries of the geometry, decompose it to the symmetry modes, separately solve the resulting SPs with reduced simulation domain size and simple PEC or PMC boundary conditions added or replacing the OP boundary conditions, and finally construct the OP solution by correctly summing the SP solutions.

We further demonstrated the symmetry-reduction method with two numerical examples considering the absorption cross-section of a single hexagonal nanowire with an incident plane wave and the absorption in a periodic nanowire array for a single incident plane wave or two interfering incident plane waves. The simulation results verified the approximately linear scaling of maximum RAM use with the size of the symmetry reduced simulation domain (the size mostly dictates the DOF in the system and hence the system matrix size) and also the reduced total solver time due to the superlinear scaling of the solver process with the reduced system matrix size. We also pointed out how, in the periodic array models, replacing the periodic Floquet boundary conditions in the OP full domain with PMC or PEC boundary conditions in SPs with the symmetry reduced domain additionally results in an even less dense system matrix and hence further reduced maximum RAM use and total solver time.

To give more perspective to these results, we can compare them with the computational cost gains in the similar symmetry reduction for FMM models of periodic arrays [18]. The FMM symmetry reduction method starts with the same concept of decomposing the OP to SPs with the symmetry modes, but then proceeds quite differently as the properties of the symmetry modes are used to find a reduced basis with which to solve the SPs (i.e., symmetry-adapted bases with fewer basis functions). The authors of those FMM studies demonstrated maximum RAM use reduction by a factor of 4 and 16 in the case of ${\mathrm{C}}_{\mathrm{s}}$ and ${\mathrm{C}}_{2\mathrm{v}}$, respectively, and computation time reduction by a factor of 4–64 depending on the symmetry and specific incidence configuration. In our case, with the models “SI no sym. red.”, “SI sym. red. 1/2” and “SI sym. red. 1/4”, the RAM use is reduced by approximately a factor of three and nine in the case of ${\mathrm{C}}_{\mathrm{s}}$ and ${\mathrm{C}}_{2\mathrm{v}}$, respectively. Moreover, if tuned properly, iterative solvers could also work for those models with no Floquet boundary conditions, which would likely reduce the maximum RAM use even further. On the other hand, even our largest reduction in solver time by a factor of approximately 10, with the “APMI sym. red. 1/4” model (with only one SP per polarization) compared to the "APMI no sym. red." model, is much more modest than the factor of 64 reported for the FMM (for normal incidence and x- or y-direction polarization). However, we would like to point out that there are several other factors to consider when determining which modeling method is the fastest and most efficient for the given problem at hand to begin with [20], and the symmetry reduction gains in FEM and FMM are hence not necessarily directly comparable. It is also worth noting that both our FEM and the FMM symmetry reduction approach also lead to the same restrictions for the allowed incidence directions in periodic models.

We would like to point out that our method is quite straightforward to implement for the case of a homogeneous background medium as well. First, it is possible to use our method as is but, if the nanostructure also has the $xy$-plane mirror symmetry in addition to the $xz$- and $yz$-plane mirror symmetries, this additional symmetry is not exploited. Second, taking all three mirror symmetry plane operations into account results in the ${\mathrm{D}}_{2\mathrm{h}}$ point group symmetry, where all the symmetry operations can be expressed as combinations of the three mirror plane symmetry operations, like is the case with the ${\mathrm{C}}_{2\mathrm{v}}$ point group symmetry and $xz$- and $yz$-plane mirror symmetry operations as discussed above. The ${\mathrm{D}}_{2\mathrm{h}}$ point group can then similarly be used to find the relevant symmetry modes and each can be assigned a different combination of PMC and PEC boundary conditions on the $xy$-, $xz$-, and $yz$-plane boundaries in the SPs with $1/8$ of the OP simulation domain. Note that the above discussed restrictions on incidence angles in the case of periodicity along the $xy$-plane would also apply.

In this work, we concentrated on the simple perpendicular mirror symmetry planes but it would be interesting to see if the other symmetries could be exploited in FEM modeling as well, like is the case with FMM [17]. The issue here is that the above FEM symmetry reduction relies on the use of PMC and PEC boundary conditions and, in a simulation domain cut with non-perpendicular mirror symmetry planes, the different combinations of PMC and PEC boundary conditions do not necessarily correspond to the symmetry modes. For example, the single hexagonal nanowire could be symmetry reduced with the ${\sigma }_{\mathrm{v}}^1$ and ${\sigma }_{\mathrm{d}}^1$ mirror symmetry planes (see Figure 1b) and the third side-boundary would be terminated with a PML. However, the boundary conditions would give only four different combinations, while there are 12 symmetry modes in the ${\mathrm{C}}_{6\mathrm{v}}$ symmetry (see ref. [25] for details). On the other hand, the symmetry reduced primitive cell in a hexagonal lattice would be bounded by three mirror symmetry planes (see Figure 2c) leading to eight different combinations of PEC and PMC boundary conditions. It could be possible that different linear combinations of the symmetry modes would correspond to the boundary conditions in these cases (instead of just one symmetry mode) and it would still be possible to find proper SPs from the OP decomposition. Therefore, this issue could be worth further research.