Features of reflection of electromagnetic waves from nanometric perforated multilayers including epsilon-near-zero metamaterials

Using the exact solutions of electromagnetic boundary problems, analytical modeling of reflection of electromagnetic waves from nanometric perforated multilayers has been carried out. New features of operation of the multilayers including the substrate or layers of epsilon-near-zero (ENZ) materials are established. Presence of the ENZ main layer or substrate leads to the quickly changing and extreme values of phase and module of amplitude reflection coefficients depending on the system parameters. The ENZ (or metallic for the thicker systems) substrate has a significant impact on the transformation of phase difference of the reflected waves. The detailed numerical analysis of the obtained results for the multilayers including silver or phase change materials (germanium antimony tellurium alloy, vanadium dioxide) components is performed. The considered reflection characteristics are reasonably “stable” to variation of the system parameters such as oblique incidence of the exciting radiation (for TE or TM polarization), possible presence of magnetic properties of the layers and effective electromagnetic anisotropy of the substrate material. The obtained results can be used to develop ultra-thin (with significantly subwavelength thicknesses) transformers of phase and amplitude of reflected radiation, holograms, metasurfaces and other nanophotonics applications.

The effective transformation of electromagnetic wave phase and amplitude at nanometric dimensions of a composite material is required for many "non-conventional" ultra-thin optical devices. For this purpose, metasurfaces realized on various physical principles are successfully used [2,4,5,8,10,14]. However, the nanofabrication of metasurfaces requires, as a rule, complex, time-consuming and expensive methods, such as focused ion beam or electron-beam lithography techniques. Moreover, it is also hard to realize large area devices (e.g. filters, holograms, metasurfaces) with these methods. Therefore, different alternative methods for both metasurfaces development and realization of their functionality for applications are actively proposed. Thus, the nanometric phase-only binary holograms (with size of 3 mm Â 3 mm and thickness of 60 nm) based on thin films of a topological insulator (antimony telluride, Sb 2 Te 3 ) on the silicon substrate have been realized and investigated in paper [16]. It was used a multilayer system "air À Sb 2 Te 3 film (including the bulk layer between two surface layers) À SiO 2 /Si substrate" in which the holes of diameter of 2 mm in the Sb 2 Te 3 film were formed with the laser ablation method. The proposed direct laser writing method for creation of the nanoholograms and experimental results [16] showed a number of possible advantages in comparison with the metasurface holograms (in the first place, more simple and technological method of production). In spite of the experimental verification of the concept of such holograms, the principle of operation of these multilayered systems required the revision and significant clarification [17].
New perspective possibilities of further decreasing the thickness of binary holograms and transformers of phase characteristics of reflected radiation were discovered using analytical modeling and optimization [17] of the perforated multilayers realized in [16] (for the case of the exciting wave normal incidence). It turns out, that the phase shift difference (PSD) of the waves reflected from the perforated multilayer (one of the main quantities, together with the reflection coefficients, determining the holographic image quality or phase transformation effectiveness) can substantially increase when using epsilon-near-zero (ENZ) MMs as the substrate of such ultra-thin layered systems. Various properties of ENZ materials are widely used in modern photonics applications [18][19][20]. In particular, ENZ substrates were used for improving the parameters of plasmonic antennas [21,22] and metafilm absorbers [23,24]. However, the PSD increase effect is non-trivial for the problems of light reflection from the perforated multilayers, and it can be useful for many (not only holographic) nanophotonics applications.
The present paper has three interconnected goals: (i) to obtain more general (in comparison with one in papers [16,17]) analytical model of the multilayer basic element (a micro-or nanohole) using the exact solutions of the corresponding electromagnetic boundary problems; (ii) to specify the features of operation of the nanometric perforated multilayers including the substrate or layers of ENZ MMs that determine new possibilities of obtaining ultra-thin devices with large PSD and reflection of electromagnetic waves; (iii) to execute a parametric optimization of the system under consideration for a number of typical and frequently used materials of the nanolayers and substrate.
In the work, new analytical solutions of the electromagnetic boundary problems for the perforated multilayers are obtained (Sect. 2). In contrast to the models [16,17], these solutions take into account oblique incidence of the exciting wave (for TE or TM polarization) on the layered system, possible magnetic properties of the layers and optical anisotropy of the substrate. The numerical modeling and detailed graphical analysis of the solutions (Sect. 3) explain the main reasons, exhibition conditions and features of the considered effect of increasing PSD of reflected waves in presence of ENZ multilayer components. The conclusion (Sect. 4) includes a brief discussion of the obtained results and summarizes the paper. The additional material excluded from the main text of the work is contained in Appendix A.

Statement of the problem and analytical model
The main investigated object in the present work is a microhole in the multilayered system under oblique incidence of an electromagnetic plane monochromatic wave (Fig. 1). The choice of such multilayer geometry is specified, in the first place, by the experimental realization and modeling of the similar system in [16,17] (at the normal incidence of the exciting wave and nonmagnetic isotropic layers). The microhole can be obtained with various methods, e.g. it can be ablated with the focused laser spot using direct laser writing method [16]. In the case of a phase-only binary hologram, the hole represents one of the hologram pixels (the number of pixels was 1,500 Â 1,500 for the holograms considered in [16]). The fabrication of the perforated multilayers as nanometric binary holograms [25][26][27] and the holographic imaging procedures (forming the source and reconstructed images) were described in detail in [16,28]. It is important that the hole can be one of the basic and simple "working elements" of not only holograms but although ultra-thin (with thicknesses much less than the wavelength) transformers of phase and amplitude of the reflected radiation [17]. In particular, it can be interesting for development of metasurfaces, perfect absorbers, nanocoatings, and other devices using ordered micro-or nanohole arrays [29,30].
When reconstructing the holographic image by the oblique incident laser beam, rays 1 and 2 ( Fig. 1a) reflect from the boundaries "input medium À Spacer 1" and "input medium À Additional layer (on the substrate)", correspondingly, interfere, and form the image. For nonholographic devices of transformation of reflected waves characteristics (e.g. for color filtering systems [6]), the image does not form, and waves 1 and 2 are used for realization of the required interference picture. The used coordinate system and geometry of the corresponding boundary problems are illustrated by Figure 1b. Axis Z is perpendicular to the layers boundaries and plane XZ is assumed to be the incidence plane. Phase multiplier exp[ik 0 (m 1 x + m 3 z) À ivt] is used in expressions for the fields where m = (m 1 , 0, m 3 ) = k/k 0 is a complex refraction vector [31], k is a wave vector (vectors are denoted with boldface symbols), k 0 = v/c is a wave number for vacuum, i 2 = À 1, r = (x, y, z) is a radius-vector. Meanwhile, the transversal component of refraction vector is supposed to be real (m 1 does not change at different boundaries in solution of the boundary problem) and the longitudinal component m 3 is complex in the general case. Layers of thicknesses d l (here and below indexes a, l = 1, 2, 3, 4, b correspond to the quantities characterizing the multilayer components in Fig. 1b) are, respectively, Spacer 1, the main layer, Spacer 2, and the additional layer (e.g. the thin SiO 2 layer on the Si substrate [16] or the substrate protective layer). Media in front of Spacer 1   Fig. 1. (a) Scheme of the multilayer with the microhole that can be obtained with the laser ablation method (one of the pixels in the case of the phase-only binary nanohologram). The rays reflected from the multilayer are denoted by numbers 1 and 2 (the incident wave is not shown). D is the hole diameter. (b) Geometry of the corresponding electromagnetic boundary problems. The multilayer includes additional (protective, spacer) "conventional" layers and the main layer and optically thick substrate which electromagnetic properties can correspond to ENZ MMs. E in1,2 and E r1,2 are the fields of the incident and reflected waves.
(e.g. air) and behind the additional layer (the substrate) are assumed to be semi-infinite (optically thick). So, the reflected wave inside the substrate is absent. As in [16,17], electric field strengths of the incident and reflected waves are considered at the top of each interface (z = 0 for ray 1, z = z 3 = d 1 + d 2 + d 3 for ray 2, Fig. 1b).
Let us investigate the systems with rather "wide" and "shallow" holes corresponding to the condition d 1 + d 2 + d 3 ≪ l < D, where l is the wavelength in the input medium, D is a certain linear size of the hole (e.g. diameter, Fig. 1a) in plane XY. As we are interested in, mainly, optical applications, the values of D are assumed to be of the order of several micrometers (D = 2 mm for the data of [16]). The same size order (several wavelengths) takes place for horizontal (in plane XY) dimensions of the structure (Fig. 1b) in the space occupied by the spacers 1, 2 and main layers (and for distances between the microholes).
Isotropic materials of the input non-absorbing medium and absorbing layers are characterized by dielectric permittivity, e a , e l , and magnetic permeability, m a , m l , correspondingly. In the general case, we assume that the substrate material is characterized by the tensors is double singular and e e (m e ) is non-singular complex eigenvalue of tensor e (m). So, the substrate is assumed to be the uniaxial magnetic medium with the optical axis which is perpendicular to the layers boundaries. Such anisotropy can take place, in particular, when the ENZ substrate is realized using layered-periodic metal-dielectric structures [32]. The scalar permittivities (permeabilities) or all or some components of tensors e b , m b can have the negative real parts, all the imaginary parts of these quantities and refraction indexes n a;l ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi e a;l m a;l p are positive according to the causality condition. So, we consider absorbing linear materials of the layers and substrate including the possible case of the MM (or metallic) layers or substrate. The denotations a 0 = Re (a), a 00 = Im (a) will be also used for the real and imaginary parts of scalar or vector quantities.
Further, we use the constitutive relations D = eE, B = mH and Maxwell's equations for monochromatic waves [31] where e À1 , m À1 are inverse to e, m tensors, and the dispersion equations (m 2 = em for isotropic media) for the materials of the layers and media outside the multilayer. Then the relations for the fields (TE and TM waves) in the isotropic media (except the substrate, Fig. 1b) can be represented in the following form (the general phase multiplier exp [ik 0 m 1 x À ivt] is omitted in the expressions for the fields below). TE waves can be determined by the expressions where lower index v = a, l points to the correspondence of a quantity to the input medium (a) or lth layer, upper indexes "+" and "À " correspond to the transmitted and reflected proper waves in an isotropic medium (these waves decay in the direction along and oppositely to axis Z, respectively). The following denotations are also used in equations (2): In a similar manner, we have for TM waves where The dispersion equations for the substrate uniaxial material take the form Using equations (4), the expressions for the transmitted proper TE and TM waves in the substrate are Þ 00 > 0). Let us use the conventional method of solution of the electromagnetic boundary problems that based on application of the continuity of the resulting fields E, H tangential components on the multilayer boundaries. Then one can obtain the following expressions for amplitude reflection coefficients of wave 1 ( Fig. 1) for TE and TM polarizations of the incident wave are the phase shifts arising on reflection of the waves corresponding to rays 1 and 2 ( Fig. 1). Equations (7) and (8) take into account exactly all the multiple reflections in the layers and represent the exact solutions of the corresponding boundary problems. In the limiting transition to the exciting wave normal incidence, non-magnetic multilayer components, and the optically isotropic substrate, equations (7) and (8) correspond to the results of papers [16,17].
The PSD arising for reflected waves 1 and 2 at the boundary z = 0 and determining the holographic image quality or phase transformation effectiveness takes the form [16,17] where d is the phase shift due to the optical path difference of wave 2 on the path z = 0, z 3 , 0 for oblique incidence of the exciting wave at incidence angle c.
Expression (9) can be also obtained when solving the boundary problem for the particular case of the multilayer in which e 1,2,3 = e a , m 1,2,3 = m a (without attraction of quantity d). With that, the wave 2 phase shift is determined at the boundary z = 0. The proposed analytical model takes into account only two main "channels" of generation of reflected waves in the system under consideration: "conventional" reflection from boundaries z = 0, z 3 . It is obvious that a part of the radiation also reflects from the holes vertical walls with the incidence angle c growth. Moreover, the model does not account the hole form, their location in plane XY, as well as possible excitation of additional electromagnetic near-field modes including the corresponding effects of extraordinary interaction of radiation with the multilayer (in presence of the ordered hole arrays under conditions of surface plasmon excitation) [33][34][35]. However, the multilayer geometry ( Fig. 1) corresponding to the experiments [16] includes the holes with deeply subwavelength "depths" (d 1 + d 2 + d 3 ≪ l) and "rather macroscopic" dimensions in plane XY (D > l). Namely that provides the possibility to use rather simple analytical models for the reflection coefficients (Eqs. (7) and (8)) without application of more difficult analytical and numerical methods of the diffraction theory (e.g. [36,37]).
We also assume (based on the good correspondence of the modeling results for the normal incidence to the experimental data [16,17]) that the obtained simplified model can be adequate under relatively small incidence angles. This assumption can be the object of an experimental verification (or a full-wave numerical analysis) in the following investigations. One of the main aims of the analysis below is to investigate the changes of the characteristics of the considered model system under oblique incidence of the exciting wave. As a result, the incidence angle regions will be determined when the multilayer still has the useful properties observed for the normal incidence. That is what shall be the main criterion of a smallness of investigated incidence angles further.

Numerical and graphical analysis
The main task of the following numerical analysis is the search of parameters and conditions of use of the perforated multilayered systems providing simultaneously large values of waves 1, 2 PSD and reflection (the main characteristics required for applications) at the minimal (significantly subwavelength, of the order of l/10 and less) multilayer thicknesses. Accounting equation (9), it means that we investigate the cases when j' TE;TM [16]). Furthermore, indexes TE, TM are omitted when both types of the waves are considered or the choice of the waves is clear from the context. Primarily, the additional possibilities of optimization of the system (in the first place, the whole thickness reduction) arising in presence of a MM (especially, ENZ) main layer or substrate are investigated. These possibilities were firstly identified in paper [17] for a number of particular cases of the considered multilayers and required an additional analysis in the general case.
Let us note that, at first glance, there is a trivial method to obtain large PSD values (of the order of p) of reflected waves 1 and 2, |' 1 À ' 2 |, from the boundary of the halfinfinity (optically thick) media. Using the expression for amplitude reflection coefficient at normal incidence, r ab = (n a À n b )/(n a + n b ), where n a,b are reflection indexes of the first (a) and the second (b) media, one can simply show that PSD of these waves is near p at reflection of wave 1 (2) from the optically more dense (less dense) medium relatively to the input medium. Moreover, this conclusion is also true in presence of some absorption in the media and under oblique incidence of the exciting wave for small incidence angles. As a result, a transformation of the boundary (that is, of medium b) into small sections where n b > n a and n b < n a provides satisfaction of the condition |' 1 À ' 2 | ≈ p even for the ideal case However, in this case, reflection under condition n b < n a is very small (|r ab | 2 < 0.1) for applications even for rather large values of dielectric contrast |n a À n b |. In addition, a violation of the requirement of optically thick medium b for reflection of any of waves 1, 2 makes the considered problem to be nontrivial and difficult for analytical solution even for the simplest case of the one main layer system (under condition d 1,3,4 = 0).
For the numerical and graphical analysis, we use some model ENZ MMs and two types of the materials frequently used in photonics applications: silver as a metallic component and dielectric phase-change materials (PCMs) that are widely used in modern optical data storage devices, optoelectronics, holography, and other applications [38][39][40][41][42][43]. PCMs undergo phase transitions resulting in drastic changes of their material (especially, dielectric) properties under changes in temperature. As a result, there are wide opportunities of active thermal and optical control of various devices on the basis of PCMs. We will consider the cases when the main layer or substrate are originated by Ag or typical PCMs such as germanium antimony tellurium alloy (Ge 2 Sb 2 Te 5 , GST 225) in amorphous and crystalline states and vanadium dioxide (VO 2 ). The literary data contain a rather wide variation of the Ag, GST 225, VO 2 permittivity values (m = 1 for these materials) subject to fabrication methods, thicknesses of the layers, etc. For the numerical analysis, the following values for silver were used: e Ag = À 4.38 + 0.20i, À12.72 + 0.52i, À23.69 + 1.07i for l = 400, 550, 700 nm, correspondingly [44]. PCMs permittivity was chosen at l = 650 nm: e GSTcryst. =9.72 + 12.96i, e GSTamorph. = 10.35 + 7.48i [41], e VO 2 ¼ 6:99 þ 1:98i [40]. For all the figures below, air is assumed to be the input medium (e a = m a = 1).
As the main layer thicknesses, we use the values d 2 = 4, 5 nm (such layers will be conditionally called "thin") and d 2 = 40, 50 nm ("thick" main layers). It is assumed for many figures below that layers Spacers 1, 2 and additional one are absent (d 1,3,4 = 0), as these layers have a rather weak impact on the investigated multilayer characteristics for thicknesses d 1,3,4 of the order of several nanometers and under satisfaction of the condition d 2 ≫ d 1,3,4 [17]. The impact of both these layers and possible anisotropy of the substrate material are discussed in the end of this section.
All the phase shifts in the graphs below are calculated in the relative units of radians/p. Phases ' 1,2 are determined below using built-up functions of the mathematical software (Mathcad, PTC) using the complex values r 1,2 calculated according to equations (7) and (8). As a rule, these functions return the principal argument of a complex number, between -p and p. As a result, phase values "leaps" can arise in graphs of dependences of quantities ' 1,2 on the system parameters at phase crossing of the real axis of the complex plane. So, for clarity, we also use the positive angle u = arccos((r 1 r 2 )/(|r 1 ||r 2 |)) (0 u p) between vectors r 1;2 ¼ ðr 0 1;2 ; r 00 1;2 Þ in the complex plane, together with ' 1,2 values.
At first, let us consider some general features of the investigated systems that take place for both oblique (under rather small incidence angles) and normal incidence of the exciting wave (when the difference between TE and TM waves disappears) independently of the presence or absence of magnetic properties of the layers. The modeling results for the choice of electromagnetic properties of the substrate near to ones of ENZ MMs and Ag, PCMs as the main layer material are given in Figures 2-6.
For both very thin and thicker main layers (the left and right panels of Figs. 2-6, correspondingly), presence of the ENZ substrate provides near to extreme values of quantities r 0;00 1;2 , ' 1,2 , u, D', |r 1,2 | 2 due to features of dependences r 0;00 1;2 ðe 0 b Þ. In particular, the extreme values of quantities ' 1 and (or) ' 2 , u, D', r 0 1;2 take place on intervals between two points of axis e 0 b corresponding to the conditions r 0 2 ¼ 0 (for both thin and thick main layers) and r 0 1 ¼ 0 (only for thin main layers) (Figs. 2 and 3a, b, d, e). For the case of the boundary of two half-infinity media (a, b) and normal incidence of the exciting wave, one can simply show that the condition r 0 2 ¼ 0 leads to the relation |n a | = |n b | for refraction indexes. This relation is satisfied twofold when parameter e 0 b varies (at the values slightly to the left and to the right of the point e 0 b ¼ 0) for the data in Figures 2 and 3. There are extremes of quantities r 0;00 1;2 on these intervals (Figs. 2 and 3a, d), whereas maxima of quantities r 0 1;2 take place for thin main layers (Figs. 2 and 3a). So, with the growth of the main layer thickness up to the values of the order of l/10 and larger, the reflection characteristics of wave 1 begin to depend weakly on the substrate parameters, and quantity r 1 is determined, substantially, by the main layer parameters. This is particularly true for the silver main layer (Fig. 2d).
The graphs of trajectories of the vectors r 1,2 tips on the complex plane (Figs. 2 and 3c, f) point to the determining role of the substrate (that is, the reflection characteristics of wave 2, Fig. 1) in realization of large PSD values under changes of parameter e 0 b . For clarity, characteristic points are chosen on these trajectories (Figs. 2 and 3c, f) that illustrate the turns of vectors r 1,2 around the origin of coordinates of the complex plane. It is seen that the tip of vector r 2 passes the much longer way (it "moves" faster with the growth of e 0 b ) in comparison with the motion of the vector r 1 tip (especially, for the thick main layers, Figs. 2 and 3f). With that, larger angles u and lengths of the vector |r 2 | are provided. The main changes of vectors r 1,2 take place within the parameters range corresponding to the ENZ substrate (Figs. 2 and 3c, f). This is particularly true for the thicker main layers originated by both Ag and PCMs (Figs. 2 and 3f).
For small values of |e b |, quantity |' 1 À ' 2 | exceeds significantly (by order of magnitude and greater, especially for thin layers) phase difference d of waves 1 and 2 corresponding to the optical paths difference of these waves (Figs. 2 and 3b, e, values of d are calculated according to Eq. (9) at n a = 1). With that, the maximal in magnitude values of angle u and PSD |D'| do not correspond to the maxima of reflection of waves 1 and 2 (Figs. 4, 5, a, f, d, e, i, k). However, under condition |e b | ≪ 1 (where is, e 0 b ≠ 0), the realization of rather large reflection of the waves and values |D'| is possible, e.g. |r 1,2 | 2 ≈ 0.5, |D'| ≈ 0.3 Ä 0.4p at d 2 = 4 nm (Fig. 4b, c, d, e), |r 1,2 | 2 ≈ 0.8 Ä 0.9, |D'| ≈ 0.8 Ä 1.0p at d 2 = 40 nm (Fig. 4g, h, i, k). In this case, the increase of the incident wave wavelength leads to the growth of PSD and reflection (Fig. 4). Let us note that in the case of the "conventional" (not ENZ) substrates for the similar systems, values of |D'| decrease with the growth of l and the condition |' 1 À ' 2 | ≪ d is satisfied, especially, with increasing the main layer thickness d 2 [16,17]. According to the data in Figures 3 and 5, the considered extreme character of the dependences of PSD and reflection for the ENZ substrate on the system parameters also takes place when PCMs are chosen as the main layer material.
Dependences of the main reflection characteristics of the multilayers on the incidence angle (under c = 0 Ä p/5 and TE, TM polarizations of the incident wave) for the parameters corresponding to Figure 4 data are illustrated by Figure 6. Here it is chosen: l = 550 nm, e b = 0.1 + 0.1i, the main layer material is silver. It is seen that the phase characteristics (D', d, |D'|/d) change rather weakly with the growth of c in the considered range. With that, the variation of reflection coefficients is stronger, especially for the thin main layers (Fig. 6a, c). The values of quantities |D' TE,TM |/d are in the ranges 7 Ä 13 (d 2 = 4 nm, Fig. 6a, c) and 2 Ä 4 (d 2 = 40 nm, Fig. 6b, d). So, the thin main layers provide a much better satisfaction of the condition |D' TE,TM | ≫ d, however the values |D' TE,TM | are significantly larger for the thicker layers (Fig. 6b, d). With the growth of c for the given nonmagnetic multilayer components, the TM-excitation provides somewhat larger values |D'| for thick layers (the data of Fig. 6d in comparison with Fig. 6b). The reflection characteristics of thin layers for TE and TM excitations are comparable (Fig. 6a, c). Thus, Figure 6 data characterize the possibilities of simultaneous realization of rather large values |D' TE,TM |, |r 1,2 | 2 for both the normal and oblique incidence of the incident radiation on the ultrathin multilayers (for rather small c values).
The similar numerical analysis also showed that the transition to the oblique incidence and TE or TM polarizations of the exciting wave changed weakly the other graphs in  to incidence angles of the order of c = p/6 Ä p/5, the graphs changes did not exceed several per cents, especially for TE waves). So, the considered effects of the PSD increase in presence of the ENZ substrate are rather "stable" to the incidence angle changes. That can be used for spatial separation of the incident and reflected waves for oblique incidence and decrease of the undesirable effects of their interference in the region in front of the multilayer.  (Figs. 7 and 8b, c). In the case of the silver substrate, rather large values of PSD and reflection, |D'| ≈ 0.3p at |r 1,2 | 2 ≈ 0.9, are possible (Fig. 7b, c). With that, we have |' 1 À ' 2 | ≫ d = 0.04p at l = 400 nm. According to the numerical estimations, the following parameters of the main layer are required in this case: e 0 2 ≈ 20 and larger, e 00 2 < 1.
For the thick main layers (the right panels of Figs. 7 and 8), the silver substrate is also more optimal in comparison with the PCMs substrate. In this case, rather large values of parameters |D'|, |r 1,2 | 2 can be obtained (Fig. 7d, e, f). These values are compatible with ones for the case of ENZ substrate (Fig. 4f, g, h).
The curves representing the pairs of values (|r 1,2 | 2 , u) corresponding to the data in Figures 4 and 5 and Figures 7  and 8 are given in Figures 9 and 10, respectively. It is seen that the maxima of values |r 1,2 | 2 , u for thin layers are not realized simultaneously (the left panels of Figs. 9 and 10). However, the usage of the silver main layer or substrate (for both thin and thick main layers, Fig. 9) provides rather large values of (|r 1,2 | 2 , u) (with that, je 0 2;b j ≪ 1, e 0 2;b ≠ 0), in contrast to the similar curves for the case of PCMs multilayer components (Fig. 10, a, c). For the thick main layers (the right panels of Figs. 9 and 10), the regions of values (|r 1,2 | 2 , u) can be possible where |r 1 | 2 , u are near to the maximal values (Fig. 9b, curves at l = 550, 700 nm,  Fig. 10b). The cases of PCMs substrates are not optimal (Fig. 10c, d). The data in Figures 9 and 10, also point to an ambiguous, in the general case, correspondence of the |r 1 | 2 and u values (when parameters e 0 2;b change) as some curves have loops.
For the data in Figures 7 and 8, the numerical analysis shows that the transition to the oblique incidence of the exciting wave (up to incidence angles of the order of c = p/6 Ä p/5) changes rather weakly the graphs only for TE polarization. Apparently, that is due to the choice of nonmagnetic multilayer components in this case. For TM polarization (the data in Fig. 11 in comparison with one in Fig. 7, the silver substrate), the additional peaks in the graphs of dependences u e 0 2 À Á , D' e 0 2 À Á (Fig. 11a, b, d, e) and the corresponding dips in graphs for jr 1 j 2 e 0 2 À Á (Fig. 11c, f) arise under the condition je 0 2 j ≪ 1. The similar features for oblique incidence take place also in Fig. 8 (for the PCMs substrate).
Properties of the considered multilayers under changes of parameter m 0 b for the silver and PCMs main layer are characterized by Figures 12 and 13  layer thickness varies. Under condition d 2 ≫ d 1,3,4 (Fig. 14a,  b), the presence of the spacer and additional layers leads to the changes of D' and |r 1 | 2 values of the order of 5-10 per cent and less. It is the case that was realized experimentally [16]. With that, one can replace the multilayer by one main layer on the substrate (that is, by more simple system) without significant changes of phase and reflection coefficient of wave 1 [17]. With the growth of d 1,3 values, the impact of the spacer layers increases and the errors can be considerable for the compatible values of thicknesses d 2 and d 1,3,4 (Fig. 14c, d). These conclusions are true for both oblique (for TM or TE waves) and normal incidence of the exciting wave, presence or absence of the layers magnetic properties.
The numerical analysis also showed that possible optical (dielectric or (and) magnetic) anisotropy of the substrate did not change (for the incidence angles up to the values of the order of c = p/6 Ä p/5), in general, the main features of reflection of the considered multilayers. The corresponding graphs (as a part of the graphs for oblique incidence of the exciting wave) are similar to ones discussed above, so they are not shown here. In particular, the values of quantities in

Conclusion
We have considered the basic features of reflection of electromagnetic waves from one microhole in ultra-thin (with significantly subwavelength thickness) perforated multilayer for which the condition |' 1 À ' 2 | ≫ d is satisfied. Namely this condition satisfaction is the main "way" to obtain the nanometric layered systems for effective transformation of the reflected radiation phase. The modeling results for "thin" main layers (their thicknesses can be estimated by the conditions d 2 10 nm and d 2 l/60 Ä l/50) characterize, apparently, the minimal possible multilayer thicknesses which can be of interest for the nanophotonics applications connected to phase transformation of reflected waves. When the  main layer thickness increases up to 40-50 nm (under conditions l/60 Ä l/50 ≪ d 2 l/10, "thick" main layers), the basic "operation conditions" of the considered systems are realized. These regimes can be realized in the presence of ENZ components as the multilayer main layer or substrate. The ENZ main layer or substrate leads to the large changes and extremal values of the multilayer reflection characteristics (phase, the real and imaginary parts of amplitude reflection coefficients of waves 1 and 2) depending on the system parameters. The ENZ substrate (that is, the wave 2 reflection characteristics) affects significantly possibilities to obtain large PSD values and rather high reflection. With that, the wave 1 amplitude reflection coefficient changes much weaker for the main layer made of various materials (topological insulator antimony telluride [16,17], silver, PCMs). That points to the possibilities  of wide choice of materials for the main and additional (spacer, protective) layers accounting requirements of applications. The investigated features of reflection from the ENZ substrate (the wave 2 characteristics, corresponding, e.g. to the data in Figs. 2 and 3) can be of independent interest for the problems of phase and amplitude transformation of the radiation reflected from various nanometric systems including metasurfaces, metaabsorbers, nanocoatings and others.
It has been ascertained that the conditions of realization of maximal values of PSD and energy reflection coefficients are not coincident for the used diapasons of parameters. However, under small "deviations" from the maximal values, the ranges of parameter values are possible when PSD and reflection of waves 1, 2 are rather large for the multilayer thicknesses up to several nanometers (e.g. for the data in Figs. 4, 5 d, e, i, k). Along with usage of MMs components of the multilayer, "conventional" metallic substrates (in the case of rather small metal absorption, e.g. for silver, Fig. 7d, e, f) can be of interest. But in this case, in contrast to the ENZ substrate, it is difficult to obtain very thin systems (one can compare the data in Fig. 7a, b, c and Fig. 4 a-k).
The considered features of the perforated multilayers are rather "stable" to variation of the parameters: for oblique incidence of the exciting wave (for TE or TM polarization, up to incidence angles of the order of c = p/6 Ä p/5 depending on the parameters values), possible presence of the layers magnetic properties and effective electromagnetic anisotropy of the substrate, the "non-idealness" of ENZ MM properties [47]. The modeling and graphical analysis results also point to the wide possibilities of dynamical control of the reflection characteristics of the considered systems. In particular, that can be realized using PCMs components as the layers or substrate (in the cases of additional thermal or optical impact). Furthermore, the obtained results can be used to develop modifications of the considered systems for various applications, e.g. using several types of the holes providing more gradations of the phase and amplitude characteristics of the reflected waves.
The author expresses the sincere appreciation to Dr Yuri Gritsai for the very fruitful discussions of this paper results. Fig. 11. The same data as in Figure 7 but for TM polarization of the incident wave and incidence angle c = p/5. Parameter e 0 2 is used as the abscissa axis for all the graphs.