Optimization of an isotropic metasurface on a substrate

. Mathematical statement of one-wavelength antireflective coating based on two-dimensional metamaterial is formulated for the first time. The constraints on geometric parameters of the structure are found. We propose a penalty function, which ensures the applicability of physical model and provides the uniqueness of the desired minimum. As an example, we consider the optimization of metasurface composed of PbTe spheres located on germanium substrate. It is shown that the accuracy of the minimization with properly chosen penalty term is the same as for the objective function without it.


Introduction
Last few years the designs of nanostructured coatings with the reflection coefficient close to zero attract a great attention. Such coatings are promising for solar cells and other photovoltaic elements which work both in the visible and in the infrared ranges. Nowadays, high refractive index all-dielectric meta-atoms are used [1], [2] instead of plasmonics [3], [4] in order to reduce Joule losses.
Commonly, the properties of substrated metasurfaces are calculated numerically. The computations are complicated due to big divergence of characteristic scales: resonator size can be 3-20 times smaller then the wavelength . Consequently, it is necessary to choose nonuniform grids with extra fine steps to describe all areas accurately. It makes computations ineffective for optimization problems. To increase the productivity, we propose to use analytical formulas from a combination of physical models [4]- [6]. However, each model has its applicability limitations. Moreover, there are restrictions on structure geometric parameters caused by fabrication limitations. They should be taken into account to obtain reasonable solutions. As a result, optimization parameters vary only in some ranges. The optimization problem should be stated as a nonlinear inverse problem of conditional minimization.
Due to resonant response of the particle array, there are numerous peaks and dips in the metasurface spectrum. Therefore, the result of objective function minimization strongly depends on the initial approximation. By performing calculations with several raffled off initial guesses, it is impossible to guarantee that the deepest of minima we found is global rather than a local one [7]. We cannot be sure that another deeper minimum does not exist. In this paper, basing on the idea of the penalty function method, we propose a well-posed statement of the inverse problem of one-wavelength antireflective coating based on isotropic two-dimensional metamaterial. The formulation allows to find global extremum, the location of which is approximately known from physical considerations. To solve the problem, we use the interior point method [8]. Its stability and accuracy are discussed.

Problem of one-wavelength antireflective coating
Depending on the specific formulation of the problem, it is required to minimize or maximize reflectance, transmittance, absorptance or their combination. Commonly, a list of materials used for fabrication of particles and a substrate is known in advance. The parameters to be determined are period of the structure and radius of the meta-atoms. According to the Sveshnikov-Ilinskiy approach [9], the solution of the optimization problem is reduced to multiple solutions of the direct problem (in our case, calculations of electrodynamic characteristics of the substrated metasurface) with directionally modified optimization parameters. To simplify calculations, it is preferable to model the structure under study by combining numerical algorithm (for the objective function minimization) with simple analytical formulas (to solve the direct problem). Similar joint approach is often used for designing multilayer coatings with the given properties [10], [11].

Physical statement of the problem
To start with, consider square periodic array composed of spherical dielectric scatterers with refractive index and radius . The one-layer structure is located at "air-dielectric" interface with refractive index of the dielectric substrate. Such isotropic metasurface (MS) with periodicity is normally illuminated by an external plane electromagnetic wave (figure 1). To describe electrodynamic properties of MS in air, A. B. Evlyukhin proposed the model of interacting induced dipoles [6]. According to the model, each sphere is replaced by a pair of electric and magnetic dipoles. To account for the interaction with other particles, Green's tensor of the medium is constructed. Such approach seems to be general since there is no homogenization of the structure [12]. In the case of normal incidence, the reflection and transmission Fresnel coefficients are where 0 = 2 / is the free-space wave number, eff e and eff m are effective electric and magnetic polarizabilities that take into account interaction between the meta-atoms in the lattice. Here and after temporal dependence is assumed to be − .
The presence of dielectric substrate influences on the field amplitude at electric (EDR) and magnetic (MDR) dipole resonances [13]. It was shown that for all-dielectric MSs, even if the refractive indexes and are high, the interaction between spherical particle and the substrate is weak enough [2]. For this reason, the MS located on the interface is modeled as imaginary sheet described with surface susceptibility electric e and magnetic m densities depending on and from (1). The reflection s and transmission s coefficients of substrated MS in the uncoupled-element model [4] are as follows: where = 2 s is a relative dielectric constant of the substrate, = 0 e / 2 and = 0 m / 2.
The above described approach gives a good qualitative description of the properties of isotropic all-dielectric MS on a substrate. Namely, it predicts the number of maxima and minima in the spectrum and dipole resonances positions [2]. This is quite enough to use it as a block for a direct problem solution. However, if more accurate model is proposed, formulas (2) will be easily replaced by the refined ones.

Constraints on geometric parameters
Limitations on structure periodicity and meta-atom size can be of several types. Firstly, there are conditions imposed from physical considerations. Obviously, the geometric parameters of the MS are positive quantities > 0 and > 0 and the particles do not touch each other > 2 . Secondly, there are limitations associated with the fabrication process. Thus, radius of identical spherical particles manufactured by dielectric material is usually not less than 50 nm. They are not located on the substrate closely, but with the interval equals to the particle diameter or more, therefore, ⩾ 4 . And, thirdly, it is necessary to take into account the conditions when the physical model works.
In our case, we should keep in mind that the Evlyukhin model gives correct results only when the dipole approximation is applicable. In [6], the condition for the minimal period is derived . ( The maximum radius for lossless materials can be found from the criterion which requires that the dipole contribution to the scattered radiation is greater than or equal to 95% [14], [15]: Some conditions listed in this subsection are overlapped. To find physical solutions, the strongest ones should be used. In addition, as upper limit on , it seems reasonable to choose ⩽ in order to exclude far-located and, thus, weakly interacting meta-atoms.

Objective functions and mathematical statement of the problem
Consider the simplest formulation of the problem: the reflectance should be minimized at some fixed wavelength = * . We introduce the vector x = { , } describing the optimization parameters. Denote the reflectivity of the structure | s (x, )| 2 as (x, ). Let 2 be a two-dimensional vector space, 2 is the closed convex set Then our goal is to determine the vector x which minimizes the function A preliminary analysis of the function (x, ) behavior shows that it strongly depends on the particle radius ( figure 2). For small , its values practically do not change. The presence of such a horizontal plateau, which is a local minimum, leads to computation looping and further breakdown. At resonant radii, there are deep "ravines". Imposing restrictions on , we exclude needless ravines: only dipole resonance is located to the left of max . However, such limitation does not eliminate the plateau. Therefore, the problem remains multi-extremal.
To make the desired minimum unique, we modify (x, ) by adding a term in the form ( ) = ( + ) , where is an even natural number. Simple estimation for the radius 0 , corresponding to MDR at the wavelength * , is known [16], [17]. The figure 3 shows a symmetric "gutter" centered at Due to ( ) selection in the form (7), we discard minima at large radii (for which the dipole approximation does not work) and make the plateau at small radii non-horizontal, see the figure 4. Thus, ( ) is a penalty function that keeps one of optimization parameters within certain range. Finally, mathematical statement of the problem is the following (8) Figure 4. Dependencies of the objective functions (8) with different values of (colored curves) and (6) (black curve) on radius for fixed period = 5.4676 m

Optimization of the structure
Calculations were carried out for substrated MS (figure 1) with = 5 (lead telluride PbTe) and = 4 (germanium Ge) at the wavelength * = 10 m, which approximately corresponds to the human body temperature. The direct problem (i.e., one of the optimization algorithm blocks) is solved using analytical formulas (1)-(2). To find global minimum of | s | 2 for ∈ [4 , 12] m and ∈ [0.05, 1.1983] m, the standard function fmincon from MATLAB Optimization Toolbox was used. It solves minimization problem of a scalar nonlinear function of multiple variables with constraints using the interior point method.

Practical recommendations
For a prevailing part of software packages, the number of function evaluations is limited by default (i.e., there is the maximum number of iterations). For example, fmincon permits only 3000 evaluations. This measure prevents cycling. However, in the case of low gradient of the objective function, it stops the calculations before some minimum is found.
In the figure 5, the percentage of the initial approximations, for which numerical calculations converge to the minimum, is indicated near the points.
For small values of , this value is 100%, and while grows it decreases. The reason for this is as follows. In the objective function (8) with the penalty term, between the steep wall for large ( − 0 ) and the minimum there are quite flat areas (minimum "sides"), which become flatter with increasing of (figure 4). These areas require more steps than available. Computations are interrupted and fmincon returns an error. In this case, it is recommended to take the last obtained values of and as new initial approximations and continue minimization. Since these sides are flat, but not horizontal, calculations converge to the minimum point.

Comparison of the objective functions
To demonstrate the advantages of our approach, we compared the results of minimization with two objective functions (6) and (8). Initial approximations were chosen randomly: 10 computations were carried out with 100 points.
Their coordinates had Gaussian distribution, the average and the standard deviation were 0 for meta-atom radius and 4 0 for structure periodicity. Each of these initial approximations was used for both objective functions.
For the minimum of | s | 2 , the dependence of the depth on is shown in the figure 5. For comparison, black line corresponding to the averaged value of minimum for the objective function (6) is added. It is clear that, using the penalty function with power ∈ [2, 6], we make the depth smaller because the center of the gutter 0 does not exactly coincide with the coordinate of minimum point ( figure 4). With the growth of , bottom of the gutter becomes flatter, and the depth increases. Starting from = 10, the minimum depth is almost independent of and does not differ from the value obtained without the penalty term (7).

Choice of the penalty function power
To choose the value of , we were guided by the following considerations. The penalty function is introduced in order to eliminate all local minima that are not located near to MDR (approximately 0 ) and EDR (close to max ) or between them. Therefore, the power should satisfy the following conditions. On the one hand, the value ( ) has to be be greater than 1 outside the specified range (all extra minima are automatically excluded from consideration). On the other hand, the penalty term should not distort the objective function (8) outside the range. These requirements are satisfied for = 10 best of all.
Note that the usage of the penalty function (7) with power = 10 practically does not affect the accuracy of obtained geometric parameters and , since it has a very flat bottom and does not distort the objective function (8). The figure 6 illustrates the accuracy of the obtained solutions using the interior point method versus the power of the penalty term ( ). The accuracy is the distance between the minimum points of the objective functions (6) and (8) under consideration = √( − ( )) 2 + ( − ( )) 2 , where and are the coordinates of the best result of minimization without the penalty function (| s | 2 ≈ 8.0579 ⋅ 10 15 ). Here ( ) and ( ) denote minimum coordinates of (8). Because of the presence of the penalty term ( ) with insufficiently at bottom at the vicinity of the desired minimum, for small values of ∈ [2,6], the accuracy, with which and are found, is not high enough (figure 6). Beginning with = 8, the accuracy of the results of minimization coincides with the tolerance of fmincon function that is 10 −6 .

Results of the minimization
The results of one of the computations with 100 random initial approximations for = 10 are depicted on the graph of | s ( , )| 2 | (figure 7). The domain of the arguments is shown by red lines. The results of minimization are marked with light dots for the objective function (8) and with dark ones for (6) without the penalty function. It is clearly seen that in the first case all 100 points converge to the same answer that is the global minimum. However, in the second case, 43 points "get stuck" on the plateau and 1 point on the horizontal area near the right boundary of the domain. They do not reach the desired minimum. To sum up, minimization of the first objective function is complicated and unstable (it depends on the choice of the initial approximation very strongly). Because of the existence of horizontal areas, in half of the cases the computations do not provide the correct answer for the position of narrow dip to be found. On the contrary, the objective function with power-law penalty term that we have constructed allows to find the desired global minimum without reference to the position of initial points. For default number of iterations, not more than 7-10 initial approximations are required.

Conclusions
The paper is devoted to the optimization of the geometric parameters of all-dielectric high refractive index isotropic metasurface placed on a semiinfinite dielectric substrate. To solve a direct problem, it is suggested to use an analytical model combining several approaches of different authors. Constraints on period of the structure and radius of spherical meta-atoms are discussed. To construct the domain of geometric parameters, technological limitations and the conditions for physical model applicability were taken into account.
For the first time, the formulation of the problem of one-wavelength antireflective substrated metasurface is proposed, based on preliminary physical considerations about the location of narrow global minimum. Using the idea of the penalty functions, we suggest new objective function, which allows to cut off all minima except the desired one: a horizontal wide region at small radii and the local minima for large particles beyond the applicability of the dipole approximation. The results of minimization with power-law penalty term and without it are compared. The choice of the power for the penalty function providing the best result of optimization is described.
The developed technique is illustrated by the example of calculating a antireflective metasurface from PbTe on a Ge substrate for a wavelength of 10 m when both materials are non-absorbent. The reflection spectrum of the structure under consideration is constructed in the range relevant for applications from 8 to 12 m. It is shown that for non-absorbing materials, zero reflection occurs between the magnetic dipole resonance and the zero reflection region of the same metasurface, but located in the air.