Particle Swarm Optimization Algorithm for Gold Nanohole Array Design †

: Gold nanohole arrays represent a class of bidimensional plasmonic metasurfaces that are suitable for realizing optical sensors. We propose a modeling approach based on a customized particle swarm optimization algorithm implemented in the commercial Ansys Lumerical FDTD software. Providing the relevant optical and morphological parameters, we obtain a set of optimized geometrical parameters tailored to the required features. Speciﬁcally, we deal with square and hexagonal arrays of air cylinders embedded in a gold layer deposited on a glass substrate engineered for surface plasmon resonance technique-based biosensors. Two structures have been deﬁned and we plan to fabricate them to validate our design routine.


Introduction
Gold nanohole arrays (GNAs) represent a class of hybrid metal/dielectric metasurfaces, made of periodically arranged air holes drilled in a thick gold film. GNAs ability to support both localized and propagating surface plasmons [1] gives them good performances for surface plasmon resonance (SPR)-oriented applications [2][3][4][5].
Several algorithms can be used to optically tune these arrays. Among them, the particle swarm optimization (PSO) algorithm [6,7] represents a powerful and versatile tool for the design of complex photonic structures [8,9] thanks to its intuitive mathematical structure, giving total control of physical phenomena. This represents a major advantage with respect to other classes of evolutionary algorithms, such as genetic algorithms [9].
In this work, we resorted to a customized PSO algorithm implemented in Ansys Lumerical FDTD [10] where the GNAs were treated as 2D gratings following the Bragg condition. The PSO was used to tune the spectral position of the minimum reflectance value at 770 nm, set as tuning wavelength, accordingly to the requirements of the project H2020 h-ALO [11]. The algorithm was then tailored to maximize the coupling of the source power into the main plasmonic resonance for two arrays: square and hexagonal. The PSO algorithm convergence and evolution were studied, and a mathematical model was developed to read its outcomes. The mismatch in the Bragg condition to the approximate estimations intrinsically carried by the FDTD method was corrected through an iterative process inspired by [12,13].

Materials and Methods
The PSO evolutionary flow is schematically described in Figure 1. As a starting point, a general GNA structure is defined in the FDTD framework. Two parameters are defined to describe the GNA geometrical features, specifically the array pitch (P) and cylinder radius (R). Two additional parameters, the velocity for both R and P, are set to manage the swarm 2 of 6 evolution. These four parameters are randomly generated by the algorithm. R and P are used to build the different 3D FDTD structures, while the velocities are used to update the R and P values, ensuring that these fall correctly within the chosen parameter space. For each (R, P) couple, defined as agent, the corresponding simulation is run and the reflectance (Refl.) and transmittance (Tran.) spectra are calculated. Finally, the algorithm evaluates the fitness function, defined as: that can be related to the coupling efficiency of the incident source in pumping the plasmonic modes. Specifically, the algorithm maximizes the FF value with the constraint of the tuning wavelength belonging to the interval (770 ± 25) nm. The allowable FF values are compared until the highest is found through all the iterations for the whole swarm, i.e., the entire collection of agents. In the end, the couple (Best R, Best P) identifies the structure with the highest and best-tuned FF. In this case "best" means the best FDTD structure that can be built by the PSO algorithm.
Mater. Proc. 2023, 14, x 2 of 6 radius (R). Two additional parameters, the velocity for both R and P, are set to manage the swarm evolution. These four parameters are randomly generated by the algorithm. R and P are used to build the different 3D FDTD structures, while the velocities are used to update the R and P values, ensuring that these fall correctly within the chosen parameter space. For each (R, P) couple, defined as agent, the corresponding simulation is run and the reflectance (Refl.) and transmittance (Tran.) spectra are calculated. Finally, the algorithm evaluates the fitness function, defined as: that can be related to the coupling efficiency of the incident source in pumping the plasmonic modes. Specifically, the algorithm maximizes the FF value with the constraint of the tuning wavelength belonging to the interval (770 ± 25) nm. The allowable FF values are compared until the highest is found through all the iterations for the whole swarm, i.e., the entire collection of agents. In the end, the couple (Best R, Best P) identifies the structure with the highest and best-tuned FF. In this case "best" means the best FDTD structure that can be built by the PSO algorithm. Process flow diagram of the particle swarm optimization algorithm (left to right). A general FDTD structure can be mathematically defined by a point in the 2D parameter space, called agent, whose coordinates are given by pitch (P) and radius (R) values. Then, in the evolution process, the algorithm seeks in the whole space the best set of parameters (Best P, Best R) to obtain an FDTD structure that maximizes the numerical value of the selected fitness function within the imposed spectral position for the tuning wavelength.
For each specific array disposition, air cylinders are constructed on a semi-infinite SiO2 substrate embedded in a gold layer. Both cylinders and the substrate are modeled as perfectly dielectric materials with a refractive index of 1 and 1.5, respectively, while the gold layer optical properties are set by the built-in Johnson and Christy data in the material database. A plane wave source, linearly polarized along the metasurface plane, is used to illuminate the structure from the substrate side at a normal rate of incidence with a spectral interval from (650 ÷ 900) nm.
Due to the symmetric properties, the boundary conditions of the FDTD box were set to be antisymmetric along the x direction, symmetric along the y direction, and a perfectly matched layer along the z direction.
After performing convergency tests, an auto non-uniform mesh type with a fivemesh accuracy was set. According to Ansys Lumerical FDTD guidelines, due to the presence of metal/dielectric interfaces, conformal variant 2 mesh refinement was selected. To improve the simulation accuracy, a mesh override with a 4 nm resolution along the x, y, and z directions was placed in correspondence with the GNA structure. Two frequency domain and power monitors, placed behind the source and above the gold layer, were used to record the Refl. and Tran. spectra in the same spectral range. Process flow diagram of the particle swarm optimization algorithm (left to right). A general FDTD structure can be mathematically defined by a point in the 2D parameter space, called agent, whose coordinates are given by pitch (P) and radius (R) values. Then, in the evolution process, the algorithm seeks in the whole space the best set of parameters (Best P, Best R) to obtain an FDTD structure that maximizes the numerical value of the selected fitness function within the imposed spectral position for the tuning wavelength.
For each specific array disposition, air cylinders are constructed on a semi-infinite SiO 2 substrate embedded in a gold layer. Both cylinders and the substrate are modeled as perfectly dielectric materials with a refractive index of 1 and 1.5, respectively, while the gold layer optical properties are set by the built-in Johnson and Christy data in the material database. A plane wave source, linearly polarized along the metasurface plane, is used to illuminate the structure from the substrate side at a normal rate of incidence with a spectral interval from (650 ÷ 900) nm.
Due to the symmetric properties, the boundary conditions of the FDTD box were set to be antisymmetric along the x direction, symmetric along the y direction, and a perfectly matched layer along the z direction.
After performing convergency tests, an auto non-uniform mesh type with a five-mesh accuracy was set. According to Ansys Lumerical FDTD guidelines, due to the presence of metal/dielectric interfaces, conformal variant 2 mesh refinement was selected. To improve the simulation accuracy, a mesh override with a 4 nm resolution along the x, y, and z directions was placed in correspondence with the GNA structure. Two frequency domain and power monitors, placed behind the source and above the gold layer, were used to record the Refl. and Tran. spectra in the same spectral range.
The simulations were run on a liquid-cooled Intel ® Core 12th generation i7-12700K (12 core) and 64 GB of DDR5-RAM and a liquid-cooled Intel ® Core 12th generation i9-12900K (16 core) with 128 GB of DDR5-RAM. An average simulation time of about 6 h was required on the i7 for the square array, while one of about 10 h was required for the hexagonal arrays on the i9. Table 1 reports the PSO parameters, together with the FF values and spectral positions, known as λ best .

Results and Discussion
In Figure 2, the agents' velocities are plotted as functions of the iteration number for the square and hexagonal arrays, panels (a) and (b), respectively, and the evolution analysis suggests that the PSO algorithm has reached convergence. This can be understood by considering the asymptotic behavior towards the point (0,0), meaning that the algorithm is no longer evolving. Panels (b) and (d) show the same concept in relation to the FF distribution in the parameter space. At convergence, the agents start swarming around the attractor within its basin of attraction, indicating that the best FF value, compatible with the constraint on the tuning wavelength, has been reached. The corresponding (Best R, Best P) set can be used to build the best FDTD structure.
Mater. Proc. 2023, 14, x 3 of 6 The simulations were run on a liquid-cooled Intel ® Core 12th generation i7-12700K (12 core) and 64 GB of DDR5-RAM and a liquid-cooled Intel ® Core 12th generation i9-12900K (16 core) with 128 GB of DDR5-RAM. An average simulation time of about 6 h was required on the i7 for the square array, while one of about 10 h was required for the hexagonal arrays on the i9. Table 1 reports the PSO parameters, together with the FF values and spectral positions, known as λbest.

Results and Discussion
In Figure 2, the agents' velocities are plotted as functions of the iteration number for the square and hexagonal arrays, panels (a) and (b), respectively, and the evolution analysis suggests that the PSO algorithm has reached convergence. This can be understood by considering the asymptotic behavior towards the point (0,0), meaning that the algorithm is no longer evolving. Panels (b) and (d) show the same concept in relation to the FF distribution in the parameter space. At convergence, the agents start swarming around the attractor within its basin of attraction, indicating that the best FF value, compatible with the constraint on the tuning wavelength, has been reached. The corresponding (Best R, Best P) set can be used to build the best FDTD structure.  (panel (a)) and the hexagonal arrays (panel (c)). In panels (b,d) we report the FF evolution in the R and P parameter space, for the square and hexagonal array, respectively. The agent number is labelled by the color.  (panel (a)) and the hexagonal arrays (panel (c)). In panels (b,d) we report the FF evolution in the R and P parameter space, for the square and hexagonal array, respectively. The agent number is labelled by the color.
At this stage, we need to go beyond the best that the PSO method can provide to precisely tune the spectral position of the minimum reflectance value at 770 nm. As a starting point, we performed a sweep on Best P. Considering a similar sweep on Best R, we obtained several (R, P) couples used to build the corresponding GNA structures and calculated the respective reflectance spectra. In Figure 3, the scatter plot of the tuning wavelengths for the different p values considered in the sweep are reported as an example for the hexagonal array.
At this stage, we need to go beyond the best that the PSO method can provide to precisely tune the spectral position of the minimum reflectance value at 770 nm. As a starting point, we performed a sweep on Best P. Considering a similar sweep on Best R, we obtained several (R, P) couples used to build the corresponding GNA structures and calculated the respective reflectance spectra. In Figure 3, the scatter plot of the tuning wavelengths for the different p values considered in the sweep are reported as an example for the hexagonal array. At this point, we performed a polynomial fit, finding the third order to be able to describe the behavior properly. It can be seen that the Best P represents the inflection point of the polynomial curve. Thus, the retuning of the FF for the best FDTD structure means that the Best P must correspond to a tuning wavelength equal to 770 nm. This indicates a shift of the whole polynomial curve. Thus, in addition, for GNAs, the main plasmonic resonance tuning depends also on R, we considered building a retuning procedure based on the ratio between R and P. From the practical point of view, a smart method of doing so in mirroring the tuning wavelength values for the P points smaller than Best P with respect to its tuning wavelength. In this way, the inflection point becomes a corner, and in its close neighborhood, the new curve can be approximated by a parabola. So, the corner point is now corresponding to the stationary point of the parabola, easier to be shifted, consequently retuning the wavelength for the best FDTD structure. For this reason, we took the p value to be 770 nm, called Ptuned, and we calculated the corresponding Rtuned as: At this point, we retuned the two arrays and the results are depicted in Figure 4, panels (a) and (d). Considering the first derivative, it was easy to identify the stationary point, shifting it in the right way by applying Equation (2) and compensating for the R dependence of the plasmonic mode by slightly adjusting the Rtuned, as shown in panels (b) and (e). The tuned R and P were used to compute the Refl., and the results of this are plotted in panels (c) and (f). It can be observed that the minimum is properly tuned at 770 nm for both the arrays. At this point, we performed a polynomial fit, finding the third order to be able to describe the behavior properly. It can be seen that the Best P represents the inflection point of the polynomial curve. Thus, the retuning of the FF for the best FDTD structure means that the Best P must correspond to a tuning wavelength equal to 770 nm. This indicates a shift of the whole polynomial curve. Thus, in addition, for GNAs, the main plasmonic resonance tuning depends also on R, we considered building a retuning procedure based on the ratio between R and P. From the practical point of view, a smart method of doing so in mirroring the tuning wavelength values for the P points smaller than Best P with respect to its tuning wavelength. In this way, the inflection point becomes a corner, and in its close neighborhood, the new curve can be approximated by a parabola. So, the corner point is now corresponding to the stationary point of the parabola, easier to be shifted, consequently retuning the wavelength for the best FDTD structure. For this reason, we took the p value to be 770 nm, called P tuned , and we calculated the corresponding R tuned as: R tuned = R best × P tuned /P best (2) At this point, we retuned the two arrays and the results are depicted in Figure 4, panels (a) and (d). Considering the first derivative, it was easy to identify the stationary point, shifting it in the right way by applying Equation (2) and compensating for the R dependence of the plasmonic mode by slightly adjusting the R tuned, as shown in panels (b) and (e). The tuned R and P were used to compute the Refl., and the results of this are plotted in panels (c) and (f). It can be observed that the minimum is properly tuned at 770 nm for both the arrays. . Results of the retuning procedure for the square (a-c) and hexagonal (d-f) arrays: evolution of the tuning wavelength (Refl. minimum) as a function of the array pitch panels (a,d)); first derivative of the parabolic curves (panels (b,e)); reflectance curves of the optimized and tuned FDTD structures (panels (c,f)).

Conclusions
A design routine for the tuning of GNAs was developed based on a custom PSO algorithm implemented in the commercial FDTD software. A mathematical model was defined to describe the evolution of the PSO and consequently retune it to reach its best solution. The design routine was applied to both square and hexagonal arrays and two corresponding well-tuned structures were obtained.   . Results of the retuning procedure for the square (a-c) and hexagonal (d-f) arrays: evolution of the tuning wavelength (Refl. minimum) as a function of the array pitch panels (a,d)); first derivative of the parabolic curves (panels (b,e)); reflectance curves of the optimized and tuned FDTD structures (panels (c,f)).

Conclusions
A design routine for the tuning of GNAs was developed based on a custom PSO algorithm implemented in the commercial FDTD software. A mathematical model was defined to describe the evolution of the PSO and consequently retune it to reach its best solution. The design routine was applied to both square and hexagonal arrays and two corresponding well-tuned structures were obtained.

Conflicts of Interest:
The authors declare no conflict of interest.