CC BY
14
Discrete & Continuous Models
#& Applied Computational Science 2022, 30 (2) 115-126
ISSN 2658-7149 (online), 2658-4670 (print) http://journals-rudn-ru/miph
Research article
UDC 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2022-30-2-115-126
Optimization of an isotropic metasurface on a substrate
Zhanna O. Dombrovskaya
Lomonosov Moscow State University, 1, bld. 2, Leninskie Gory, Moscow, 119991, Russian Federation
(received: March 12, 2022; revised: April 18, 2022; accepted: April 19, 2022)
Abstract. 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.
Key words and phrases: antireflective coating optimization, penalty function method, constrains on geometric parameters, all-dielectric metasurface on a substrate
1. 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 A. 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.
© Dombrovskaya Z.O., 2022
This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/
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.
2. 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 p of the structure and radius r 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 meta-surface) 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].
2.1. Physical statement of the problem
To start with, consider square periodic array composed of spherical dielectric scatterers with refractive index n and radius r. The one-layer structure is located at "air-dielectric" interface with refractive index ns of the dielectric substrate. Such isotropic metasurface (MS) with periodicity p is normally illuminated by an external plane electromagnetic wave (figure 1).
Figure 1. Schematic representation of a metasurface consisting of spherical particles on a semi-infinite substrate. The structure is normally irradiated by an electromagnetic wave
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 R and transmission T Fresnel coefficients are
r> _ / eff „,eff\ rp _ , / ^eff I ~,eff\ /1 \
K=2p2 -am), 1 ~1+2p2 '"m^ (1)
where k0 = 2n/X is the free-space wave number, ff and a^f 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 e-zwt.
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 n and ns 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 xe and magnetic xm densities depending on R and T from (1). The reflection Rs and transmission Ts coefficients of substrated MS in the uncoupled-element model [4] are as follows:
(1 + e)(1 — ,/em) — (y/e — e) (1 + m)
Rs =
t =
(1 — e) (1 — y/em) + (e — e) (1 — m) (1 + e)(1 + m) + (1-e)(1 + m)
(2)
(1 — e) (1 — y/em) + (e — e) (1 — m)'
where e = n:: is a relative dielectric constant of the substrate, e = ik0xe/ 2 and m = ik0xm/ 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.
2.2. 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 p > 0 and r > 0 and the particles do not touch each other p > 2r. 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, p > 4r. 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
Ueff|2 + Ueff |2
n = I |ae 1 T |am 1 (3)
Pmin Im (of) + Im (aff)' (3)
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]:
r - A (4)
Some conditions listed in this subsection are overlapped. To find physical solutions, the strongest ones should be used. In addition, as upper limit on p, it seems reasonable to choose p < A in order to exclude far-located and, thus, weakly interacting meta-atoms.
2.3. Objective functions and mathematical statement of the problem
Consider the simplest formulation of the problem: the reflectance should be minimized at some fixed wavelength A = A^. We introduce the vector x = {p, r] describing the optimization parameters. Denote the reflectivity of the structure lRs(x, A)|2 as f(x, A). Let E2 be a two-dimensional vector space, C2 is the closed convex set
C2 = {x £ E2 ■ 4r < P < A max fbric < r < ^]- (5)
Then our goal is to determine the vector x which minimizes the function
f(x,\) = min, x £ C2• (6)
A preliminary analysis of the function f(x, A) behavior shows that it strongly depends on the particle radius (figure 2). For small r, 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 r, we exclude needless ravines: only dipole resonance is located to the left of rmax. However, such limitation does not eliminate the plateau. Therefore, the problem remains multi-extremal.
To make the desired minimum unique, we modify f(x, A) by adding a term in the form y(r) = (Ar + B)@, where P is an even natural number. Simple estimation for the radius r0, corresponding to MDR at the wavelength A^, is known [16], [17]. The figure 3 shows a symmetric "gutter" centered at r0 — X^/(2n) with width equal to (rmax — r0). Changing the value of 0, it is possible to control the slope of walls and the flatness of its bottom.
12 10 8
0.8
0.6
4
2 0
ñ
0.4
0.2
0 0.5 1 1.5 2
r, ¿¿m
Figure 2. Dependence of |Rs|2 on period p and radius r at A = 10 ^m
Figure 3. Penalty function y(r) for different values of power (3
Assume that y(r0) = 0 and y(rmax) = 1 on the wall of the gutter, then the proposed term is
y(r) =
r — r0
^max
(7)
Due to y(r) 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, y(r) is a penalty function
that keeps one of optimization parameters within certain range. Finally, mathematical statement of the problem is the following
F?[x] = /(x,A) + y(r), x GC2. (8)
r, ¡im
Figure 4. Dependencies of the objective functions (8) with different values of (3 (colored curves) and (6) (black curve) on radius r for fixed period p = 5.4676 ^m
3. Optimization of the structure
Calculations were carried out for substrated MS (figure 1) with n = 5 (lead telluride PbTe) and ns =4 (germanium Ge) at the wavelength A^ = 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 p £ [4r, 12] ^m and r £ [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.
3.1. 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 ft, this value is 100%, and while ft 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 (r — r0) and the minimum there are quite flat areas (minimum "sides"), which become flatter with increasing of ft (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 p and r as new initial approximations and continue minimization. Since these sides are flat, but not horizontal, calculations converge to the minimum point.
-3 «100%
>100%
Dd
F -9
CD
-12
56%
-15
2 4
>99%
97%*^
\ 90% 86% t
94%""' ' 86%
6 8 10 12 14 16
ft
Figure 5. Minimal values of |Ks|2 at logarithmic scale for computations with proposed penalty function (dotted curve) and without it (straight line). The percents of initial approximations, for which fmincon function converges to the minimum (see text), are
indicated nearby
3.2. 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 r0 for meta-atom radius and 4r0 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 ft 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 ft £ [2,6], we make the depth smaller because the center of the gutter r0 does not exactly coincide with the coordinate of minimum point (figure 4). With the growth of ft, bottom of the gutter becomes flatter, and the depth increases. Starting from ft = 10, the minimum
depth is almost independent of ft and does not differ from the value obtained without the penalty term (7).
3.3. Choice of the penalty function power
To choose the value of ft, 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 r0) and EDR (close to rmax) or between them. Therefore, the power ft should satisfy the following conditions. On the one hand, the value y(r) 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 ft = 10 best of all.
Note that the usage of the penalty function (7) with power ft = 10 practically does not affect the accuracy of obtained geometric parameters p and r, since it has a very flat bottom and does not distort the objective function (8). The figure 6 illustrates the accuracy 5 of the obtained solutions using the interior point method versus the power ft of the penalty term y(r). The accuracy is the distance between the minimum points of the objective functions (6) and (8) under consideration 5 = — p(ft))2 + (r — r(ft))2, where p and r are the coordinates of the best result of minimization without the penalty function (lRs|2 « 8.0579 • 1015). Here p(ft) and r(ft) denote minimum coordinates of (8).
0
Figure 6. Accuracy of the minimization for different values of the penalty function power
Because of the presence of the penalty term y(r) with insufficiently at bottom at the vicinity of the desired minimum, for small values of ft £ [2,6], the accuracy, with which p and r are found, is not high enough (figure 6). Beginning with ft = 8, the accuracy of the results of minimization S coincides with the tolerance of fmincon function that is 10-6.
3.4. Results of the minimization
The results of one of the computations with 100 random initial approximations for fl = 10 are depicted on the graph of lRs(p,r)|21 (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.
0.2 0.4 0.6 0.8 1 1.2
r, ¿¿m
Figure 7. Results of the minimization of substrated metasurface at the wavelength A = 10 ^m. Found minima of the objective functions (8) and (6) are pointed out on the graph of the reflectance lRs(p, r)|2 with white and black markers, respectively. The number of points is indicated beside them. Red lines are the boundaries of C2
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.
4. 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 antire-flective 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 an-tireflective 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.
References
[1] K. V. Baryshnikova, M. I. Petrov, V. E. Babicheva, and P. A. Belov, "Plasmonic and silicon nanoparticle anti-reflective coatings", Scientific Reports, vol. 6, p. 22136, 2016. DOI: 10.1038/srep22136.
[2] V. E. Babicheva, M. I. Petrov, K. V. Baryshnikova, and P. A. Belov, "Reflection compensation mediated by electric and magnetic resonances of all-dielectric metasurfaces", Journal of the Optical Society of America B, vol. 34, no. 7, pp. D18-D28, 2017. DOI: 10.1364/J0SAB.34.000D18.
[3] H. A. Atwater and A. Polman, "Plasmonics for improved photovoltaic devices", Nature Materials, vol. 9, pp. 205-213, 2010. DOI: 10.1038/ nmat2629.
[4] M. Albooyeh, D. Morits, and C. R. Simovski, "Electromagnetic characterization of substrated metasurfaces", Metamaterials, vol. 5, pp. 178205, 2011. DOI: 10.1016/j.metmat.2011.08.002.
[5] Z. O. Dombrovskaya et al., "Inverse problem for recovering of meta-atom characteristics by transmittance and reflectance of a metafilm", Bulletin of the Russian Academy of Sciences: Physics, vol. 79, no. 12, pp. 14961498, 2015. DOI: 10.3103/S1062873815120151.
[6] A. B. Evlyukhin et al., "Optical response features of Si-nanoparticle arrays", Physical Review B, vol. 82, p. 045 404, 2010. DOI: 10. 1103/ PhysRevB.82.045404.
[7] N. N. Kalitkin and E. A. Al'shina, Numerical Methods [Chislennye metody], book 1, in Russian. Moscow: Akademiya, 2013.
[8] R. H. Byrd, M. E. Hribar, and J. Nocedal, "An interior point algorithm for large-scale nonlinear programming", SIAM Journal on Optimization, vol. 9, no. 4, pp. 877-900, 1999. DOI: 10.1137/S1052623497325107.
[9] A. G. Sveshnikov and A. S. Ilinskiy, "Design problems in electrodynamics [Zadachi proyektirovaniya v elektrodinamike]", in Russian, Proceedings of the USSR Academy of Sciences, vol. 204, pp. 1077-1080, 1972.
[10] V. B. Glasko, A. N. Tikhonov, and A. V. Tikhonravov, "On the synthesis of multilayer coatings [O sinteze mnogosloynykh pokrytiy]", USSR Computational Mathematics and Mathematical Physics, vol. 14, p. 135, 1974, in Russian.
[11] A. V. Tikhonravov et al., "Design and production of antireflection coating for the 8-10 ^m spectral region", Optics Express, vol. 22, pp. 32174-32179, 2014. DOI: 10.1364/0E.22.032174.
[12] Z. O. Dombrovskaya and A. V. Zhuravlev, "Investigation of the possibility of metafilm modeling as a conventional thin film", Applied Physics A, vol. 123, p. 27, 2017. DOI: 10.1007/s00339-016-0642-2.
[13] A. E. Miroshnichenko et al., "Substrate-induced resonant magneto-electric effects with dielectric nanoparticles", ACS Photonics, vol. 2, pp. 1423-1428, 2015. DOI: 10.1021/acsphotonics.5b00117.
[14] G. V. Belokopytov and A. V. Zhuravlev, "Dipole polarizability of spherical particles [Dipol'naya polyarizuyemost' sfericheskikh chastits]", in Russian, Physics of Wave Processes and Radio Systems, vol. 2, pp. 4149, 2008.
[15] Z. O. Dombrovskaya et al., "Phonon-polariton meta-atoms for far infrared range", Physics of Wave Phenomena, vol. 24, pp. 96-102, 2016. DOI: 10.3103/S1541308X16020023.
[16] A. I. Kuznetsov et al., "Magnetic light", Scientific Reports, vol. 2, p. 492, 2012. DOI: 10.1038/srep00492.
[17] D. G. Baranov et al., "All-dielectric nanophotonics: the quest for better materials and fabrication techniques", Optica, vol. 4, no. 7, pp. 814-825, 2017. DOI: 10.1364/0PTICA.4.000814.
For citation:
Z. O. Dombrovskaya, Optimization of an isotropic metasurface on a substrate, Discrete and Continuous Models and Applied Computational Science 30 (2) (2022) 115-126. DOI: 10.22363/2658-4670-2022-30-2-115-126.
Information about the authors:
Dombrovskaya, Zhanna O. — Candidate of Physical and Mathematical Sciences, Junior Researcher of Faculty of Physics, Lomonosov Moscow State University (e-mail: dombrovskaya@physics.msu.ru, phone: +7(495)9393310, ORCID: https://orcid.org/0000-0003-0609-1065)
УДК 519.872:519.217
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf
DOI: 10.22363/2658-4670-2022-30-2-115-126
Оптимизация изотропной метаповерхности
на подложке
Ж!. О. Домбровская
Московский государственный университет им. М. В. Ломоносова, Ленинские горы, д. 1, стр. 2, Москва, 119991, Россия
Аннотация. Впервые приведена математическая формулировка одноволнового безотражательного покрытия на основе двумерного метаматериала. Найдены ограничения на геометрические параметры конструкции. Предложена штрафная функция, которая обеспечивает применимость физической модели и обеспечивает единственность искомого минимума. В качестве примера рассмотрена оптимизация метаповерхности, состоящей из сфер РЬТе, расположенных на германиевой подложке. Показано, что точность минимизации с правильно выбранным штрафным термином такая же, как и для целевой функции без него.
Ключевые слова: оптимизация безотражательного покрытия, метод штрафной функции, ограничения на геометрические параметры, диэлектрическая метаплёнка на подложке
