## Abstract

The analysis and synthesis of metasurfaces are important because of their emerging applications in a broad range of the operational wavelengths from microwaves to the visible light spectrum. Moreover, in many applications, like optical nanoantennas, absorbers, solar cells, and sensing, the presence of a substrate is apparent. Therefore, understanding the effects of substrates upon the metasurfaces is important, as the substrates typically affect the resonance behaviors of particles, as well as the interactions between them. In order to consider the impacts of substrates, this paper develops a method for the characterization and homogenization of substrated metasurfaces. This approach is based on independent studies of the electromagnetic behavior of the constituting nanoparticles, and the interactions between them. It uses image theory to calculate the interaction constant tensors in the presence of a dielectric substrate. Then, the contributions of the quasi-static interaction fields of the primary and image dipoles are considered as a homogeneous sheet of surface polarization currents. Finally, the closed-form expressions for the interaction constant tensors are derived. To show the accuracy of our proposed approach, the numerical results of the method are compared to other approaches, as well as with those generated by a commercial EM solver, which are all found to be in good agreement. Moreover, the effects of the refractive index of the substrate, the geometric characteristics of the particle, and periodicity of the array are also investigated on the interaction constants. We believe that this methodology is general and useful in the design and analysis of substrated metasurfaces for various applications.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The problem of light scattering from planar periodic structures has a long history. In the literature, different solutions have been presented to study this problem [1–9]. However, they have been applied with various limitations and assumptions that, on some occasions, have resulted in findings that are substantially different from experimental or full-wave simulator data [10–12]. Since introducing the concept of metasurfaces by Kuester and Holloway, more attention has been concentrated on studying these structures [13–16]. Metasurfaces are surface versions of metamaterials. More precisely, metasurfaces are surface distributions of appropriately selected small resonant scatterers, which are becoming a serious competitor for 3D metamaterials.

The analysis of a metasurface is usually referred to as the calculation of its equivalent parameters [11,12,17]. The characteristic parameters of a metasurface include its surface susceptibility tensors, collective polarizability tensors, and impedances. The importance of the electromagnetic characterization of metasurfaces is that, once obtained, they can be considered as an equivalent of the metasurface in other electromagnetic problems.

Various approaches have been proposed to analyze metasurfaces; where have mostly concentrated on metasurfaces located in homogeneous media [18,19]. In many applications, the presence of a substrate is inevitable [20–25], as it can improve the structure’s mechanical robustness, and makes it more favorable for practical applications. Furthermore, the presence of a substrate strongly affects the resonance behaviors of particles, as well as the interactions between their electric and magnetic resonances. Therefore, analyzing metasurface in the presence of a substrates is worth investigating.

In [11], the characteristic parameters of some substrated metasurfaces have been retrieved from their reflection and transmission coefficients. The reflection / transmission coefficients of a metasurface are defined in terms of its complex susceptibilities through Generalized Sheet Transition Condition (GSTC). Then, by solving the reverse problem, the complex susceptibilities are obtained from the measured or simulated reflection / transmission coefficients. Due to the complexity of solving the reverse problem, this method has only really found application in some special cases. Moreover, the retrieval of parameters from this method does not provide independent information about the inclusions themselves, or their arrangements. Furthermore, if either the array period or the inclusion characteristics are changed, the parameter retrieval should be repeated, separately, for the new case.

Hence, an alternative approach was proposed to calculate the characteristic parameters of a metasurface based on an independent study of the electromagnetic behaviors of the constituting nanoparticles of a metasurface and the interactions between them [26–28]. However, this method has, to date, been restricted to metasurfaces in homogenous media; so, in this paper, this approach has been extended to derive the characteristic parameters of metasurfaces located at the interface of two semi-infinite dielectric media. This has been done based on the notion that the presence of a substrate modifies the properties of the resonant particles, as well as the interactions between them.

In this approach, the electromagnetic behavior of a nanoparticle can be described by its individual polarizability tensors. The individual polarizability tensors link the induced dipole moments to the local electromagnetic fields. In the presence of a substrate, the effects of the scattered fields must be added to the incident fields when calculating the local fields [29]. Moreover, in the calculation of the local fields for a nanoparticle in close proximity to other nanoparticles and in the presence of a substrate, image theory is used. In this regard, we take into account the contribution of the quasi-static interaction fields of other primary dipoles, and also image dipoles as a homogeneous sheet of surface polarization currents.

This paper is structured as follows: In Section 2, we explain the proposed approach for the homogenization and characterization of metasurfaces located on a semi-infinite dielectric substrate. Then, this method is extended to the case where a metasurface located on a finite-thickness substrate, or on a metallic ground plate separated by a dielectric spacer, using the multiple-reflection method. After presenting the analytical framework, in Section 3, some different structures are considered to verify the homogenization process — beginning with a calculation of the effective polarizability tensors of a 2D periodic array of silicon nano-spheres located in free space. After this initial verification, the effects of different parameters on the interaction constant tensors are considered. Next, a 2D periodic array of nano-sphere particles is considered on a semi-infinite dielectric substrate; and the effective polarizabilities and reflection / transmission coefficients extracted from the proposed method are presented and compared to a full wave simulation. Finally, the paper is concluded, in Section 4.

## 2. Theory

#### 2.1 Metasurface on a half-space substrate

Let us begin by considering a general metasurface sandwiched between two dielectric media with permittivities, ${\varepsilon _1}$ and ${\varepsilon _2}$ (see Fig. 1), which is illuminated by a normally incident plane wave. In general, this incident field induces both the electric and magnetic dipole moments on each nanoparticle in the array. It is worth clarifying that, in this paper, whenever referring to a metasurface, we mean a periodic planar array with similar electrically small nanoparticles and the induced dipole moments are calculated at the geometrical center of the particles. In such cases, the only significant moments are the electric and magnetic dipoles, and any other higher-order multipoles can be neglected. This condition is necessary to use the method presented in this paper. To analyze the electromagnetic responses of a metasurface to an external excitation, the scattered fields from the metasurface must be calculated. These fields can be obtained using induced electric and magnetic surface current densities, which are dependent on both the electromagnetic behavior of the constituting nanoparticles, and the interactions between them.

In fact, as shown in Fig. 2, the analysis of a metasurface can be divided into two sections: considering the electromagnetic behavior of the constituting nanoparticles by calculating the individual polarizabilities of the nanoparticle; and calculating the interaction constants, which considers the arrangement of the array. By combining these two sections, the collective polarizability tensors of a metasurface can be obtained, which is the characteristic parameters of the metasurface. In the following, we study each of these sections in metasurface analysis.

### 2.1.1 Calculation of individual polarizability tensors of nanoparticles

As mentioned before, since the constituent particles of the MSs are smaller than the wavelengths of the applied fields in the surrounding media, and the induced moments are calculated at the geometrical center of the particles, the electric and magnetic dipolar moments are deemed to be sufficient to describe the resultant scattered fields. Therefore, by definition, the induced dipole moments (${\textbf p}$ and ${\textbf m}$) can be expressed in terms of individual polarizabilities of an arbitrary bianisotropic nanoparticle as follows [26]:

Where, ${\bar{\bar{\alpha }}^{ee}},\begin{array}{{c}} {} \end{array}{\bar{\bar{\alpha }}^{mm}},\begin{array}{{c}} {} \end{array}{\bar{\bar{\alpha }}^{em}}$, and ${\bar{\bar{\alpha }}^{me}}$ are the electric, magnetic, magneto-electric, and electro-magnetic polarizabilities, respectively. The polarizability tensors are functions of the physical characteristics, material, and surrounding environment of the nanoparticle. ${{\textbf E}^{loc}}$ and ${{\textbf H}^{loc}}$ are the local electric and magnetic fields acting on the nanoparticle, respectively. When a nanoparticle is placed in a homogeneous medium, the local fields are equal to the incident fields. However, in the presence of a substrate, it is necessary to consider its effects on the nanoparticle, when calculating the local fields.

Various analytical and numerical methods have been proposed to calculate the polarizability tensors of nanoparticles [17,30–35]. The analytical methods are typically limited to simple configurations and their complexities grow considerably when applied to many engineered structures [26,35,36]. Therefore, semi-analytical and numerical techniques are appropriate for considering the more complex structures. But, most of these methods compute the polarizability tensors of a nanoparticle in a homogeneous medium. However, when the nanoparticle is located at the boundary between two different media, the problem becomes more complex. Indeed, the effect of the environments must be taken into account, when calculating the polarizability tensors. Recently, we have proposed a method for calculating the polarizability tensors of an arbitrary particle on a semi-infinite dielectric substrate [29]. In [29] the effect of the scattered fields is added to the incident fields when calculating the local fields.

### 2.1.2 Calculation of interaction coefficients

As mentioned before, when a nanoparticle is considered independently, the local field acting on the nanoparticle is equal to the incident field in the free space, and superposition of the incident field and the scattered field in the presence of the substrate. However, when considering one nanoparticle next to other nanoparticles, the effects of the other nanoparticles must be considered as well, which can be done using interaction constants [26]. In this paper, we extend the approach used in [26] to the case of metasurfaces placed in the presence of semi-infinite dielectric substrates. To do so, we use the image theorem to model the effects of the substrates (see Fig. 3).

Considering similar constituent nanoparticles in a metasurface, the local fields acting on the nanoparticles would be identical. Hence, let us consider an arbitrary nanoparticle as a reference, and position the Coordinate system’s origin at its center. In this case, the local fields acting on the reference nanoparticle would be the superposition of the incident fields (${{\textbf E}^{inc}}$ and ${{\bf H}^{inc}}$) and the fields created by the other nanoparticles at the location of the reference nanoparticle. This can therefore be written as follows:

Where, ${{\textbf E}^{int,primary}}$ (${{\bf H}^{int,primary}}$) represents the interaction electric (magnetic) field created by the primary dipoles of the array, and ${{\textbf E}^{int,image}}$ (${{\bf H}^{int,image}}$) represents the interaction electric (magnetic) field created by the image dipoles that correspond to the effects of the substrate. Hence, for a metasurface located in a homogeneous medium, ${{\textbf E}^{int,image}}$ (${{\bf H}^{int,image}}$) would equal zero. To calculate the interaction fields, since all nanoparticles assumed to be similar and excited by a normally incident plane wave, the induced dipole moments on all the nanoparticles would be identical. Therefore, by summing up all the radiation fields due to the primary dipoles created at the location of the reference nanoparticle, ${{\textbf E}^{int,primary}}$ could be calculated. This is done by the method performed in [26], which considered a circle of radius, *R*, around the reference dipole, and replaced the discrete dipoles’ distribution outside the circle with a continuous distribution of dipoles. The result is as follows:

Where $\bar{\bar{\beta }}_0^{ee}$ and $\bar{\bar{\beta }}_0^{mm}$ are free space electric and magnetic interaction constant tensors, which may be expressed as [26]:

To calculate the second part of the Eq. (3), according to the image theorem, the image dipole moments are considered to have a magnitude of $\frac{{{\varepsilon _2} - {\varepsilon _1}}}{{{\varepsilon _2} + {\varepsilon _1}}}$- times the primary dipole, and depending on the direction of the primary dipole, the image dipole could be in phase or out of phase. We can now define a parameter, named contrast factor, which is equal to the ratio of the image dipole magnitude to the primary dipole magnitude, as ${C_i} = \frac{{{\varepsilon _2} - {\varepsilon _1}}}{{{\varepsilon _2} + {\varepsilon _1}}}$, and considering this contrast factor, ${{\textbf E}^{int,image}}$ can then be calculated. Similar to the calculation of ${{\textbf E}^{int,primary}}$, we consider a circle of radius, *R*, around the image of reference nanoparticle, and replace the discrete image dipoles distribution outside the circle through a continuous distribution of the dipoles. Next, the contribution of the image of reference dipole must be added to the interaction fields. Therefore, we have:

Where ${{\textbf p^{\prime}}_0}$ and ${{\bf m^{\prime}}_0}$ are the image of the reference electric and magnetic dipoles, respectively. The vector, ${\bf h}$ is the distance between the primary dipole and its image, which is parallel to the z-axis. To clarify the notation here, we can illustrate with an example. Note that ${{\textbf E}^{sheet}}({\textbf p^{\prime}},{\textbf r^{\prime}})$ is the electric field created by the continuous distribution of the image electric dipoles, with a hole of radius, *R*.

Where ${s_0}$ is the unit cell area. Moreover, ${{\textbf E}_{p^{\prime}}}({\textbf r^{\prime}})$ is the re-radiated electric field for image of the induced electric dipoles that is expressed as follows:

Where ${\varepsilon _0}$ and ${\textbf n}$ are the free space permittivity and the unit vector along the observation vector, respectively. To calculate the integral in Eq. (7), the electric dipole moments are divided into components which are parallel (${\parallel}$) and orthogonal ($\bot$) to the array planes. The integral has also been evaluated in [37], which are reduced to simple combinations of elementary functions as:

Where ${{\textbf p^{\prime}}_\parallel }$ and ${{\textbf p^{\prime}}_ \bot }$ are image of electric dipoles parallel and orthogonal to the array plane, respectively. According to the image theorem, we have ${{\textbf p^{\prime}}_\parallel } ={-} {C_i}\begin{array}{{c}} {{{\textbf p}_\parallel }} \end{array}$ and ${{\textbf p^{\prime}}_ \bot } = {C_i}\begin{array}{{c}} {{{\textbf p}_ \bot }} \end{array}$. Moreover, ${\textbf E}({{\textbf p^{\prime}}_0},{\textbf h})$ in Eq. (6) is the electric field created by the image of the reference electric dipole at a distance, *h*, from its center. This field is calculated by substituting ${\textbf n} = {\textbf z}$ and $r^{\prime} = h$ in Eq. (8). The result is:

The calculations of other terms in Eq. (6) and also interaction magnetic field created by the image dipoles (${{\bf H}^{int,image}}$) can similarly be performed (See Appendix A). By calculating ${{\textbf E}^{int,image}}$ and ${{\bf H}^{int,image}}$, these interaction fields can be expressed in terms of the dipole moments of the reference nanoparticle, as follows:

Where $\bar{\bar{\beta }}{^{\prime ^{ee}}}$, $\bar{\bar{\beta }}{^{\prime ^{mm}}}$, $\bar{\bar{\beta }}{^{\prime ^{em}}}$, and $\bar{\bar{\beta }}{^{\prime ^{me}}}$ are the electric, magnetic, magneto-electric, and electromagnetic image interaction constant tensors, respectively; which are defined as follows:

In which *a* and $\omega$ are the array period and angular frequency, respectively. Moreover, for a square lattice, as consider in this paper, we have ${R_{}} = \frac{a}{{1.348}}$ [26]. The approximate analytical formula which uses the static solution to determine the equivalent hole radius *R* gives good results up to $ka \approx 3$[26]. Notice that, in contrast to the calculation in free space (as see in [26] and also Eq. (4)), the electromagnetic and magneto-electric interaction constants, here, are non-zero due to the effects of the substrate. Finally, by substituting Eq. (4) and Eq. (11) in Eq. (3), the total interaction constants are:

Where $\beta _{ii}^{jj} = \beta _{0ii}^{jj} + \beta _{ii}^{\prime^{jj}}$, *i *= *x*, *y*, *z* and *j *= e or m.

In order to achieve an initial validation of the proposed formulation of the interaction constants, we consider the example of a special case. In this case, we assume that a metasurface is located on a substrate with ${\varepsilon _r} = 1$, which leads to a zero contrast factor. Therefore, all of the $\beta ^{\prime}$ in Eq. (12) would be zero. As a result, the interaction constants for this special case, would be equal to the same equations in [26].

### 2.1.3 Calculation of reflection and transmission coefficients

Up to now, the method of calculating individual polarizability tensors of nanoparticles and calculating the interaction constants are explained in sub-sections 2.1.1 and 2.1.2, respectively. At this point, everything is ready for the analysis of the metasurface located on a semi-infinite dielectric substrate. By knowing the individual polarizabilities and the interaction constants, the collective polarizabilities of the substrated metasurface can be calculated. In this regard, Eq. (4) and Eq. (11) are first replaced in Eq. (3). Then, the results of this as well as Eq. (2) are substituted into Eq. (1) and averaging on the unit cell area gives the surface polarization densities in terms of the incident fields:

Where $\bar{\bar{\hat{\alpha }}}$ is the collective polarizability tensor of the substrated metasurface, and calculated as follows:

Where $\bar{\bar{I}}$ is idem factor. Briefly, in the proposed approach in this paper, for calculating the characteristic parameters of metasurfaces placed on semi-infinite substrates, the individual polarizabilities are firstly calculated and transformed to the collective polarizabilities using the interaction constants. Then, by knowing the collective polarizabilities [from Eq. (15)] which are related to the induced electric and magnetic surface current densities, one can obtain the scattered fields from the metasurface, as follows:

#### 2.2 Metasurfaces on finite thickness substrates

Up to now, metasurfaces located on an infinite dielectric substrate have been analyzed. Now, using the results obtained so far, we can analyze more practical cases, such as when the metasurface is located on a finite thickness substrate (Fig. 4(a)) or on a metallic ground plate separated by a dielectric spacer (Fig. 4(b)). To do so, the multiple reflection method can be used [38]. According to Fig. 4(a), we have the following:

In [20,39], these equations were used to analyze a terahertz absorber. However, they calculated the reflection and transmission coefficients (${r_{12}}$, ${t_{21}}$, ${r_{21}}$, and ${t_{12}}$) where the metasurface was assumed to be sandwiched between two semi-infinite media, using a full-wave simulation. However, according to the current paper, one could calculate the reflection and transmission coefficients of the metasurface in the presence of a semi-infinite dielectric substrate. The advantage of using our proposed method, over full-wave simulators, is that by changing the periodicity of the metasurface there would be no need to repeat the whole analysis process. Indeed, in our proposed approach, once the polarizability tensors of the individual nanoparticles have been obtained, there is no need to perform any full-wave simulation in the calculation of the scattered fields.

## 3. Results and discussion

In the previous section, we described our method for the homogenization and characterization of a metasurface, either in a homogeneous medium, or on a dielectric substrate. To verify the accuracy and validity of the homogenization process, this section compares the results of the study to other approaches.

Firstly, we investigate a special case in which a 2D periodic array of silicon nano-spheres is located on a substrate with ${\varepsilon _r} = 1$. As mentioned before, in the proposed homogenization method, we first need to obtain the polarizability tensors of an individual nanoparticle. In this case, we use the analytical Mie theory to calculate the polarizability tensors of the sphere, located in free space [40]. Then, by calculating the interaction coefficients, the characteristic parameters of the metasurface are obtained according to Eq. (15). Figure 5(a) shows the result of the effective polarizabilities of this metasurface. Moreover, the reflection and transmission coefficients of the proposed approach are compared to those obtained from the analytical theory presented in [41], and the full-wave simulation [42]. As can be seen in Fig. 5, our proposed approach is well in agreement with the analytical method presented in [41].

After this initial verification, we now, investigate the effects of different parameters on the interaction coefficients. Firstly, let us consider an array of arbitrary nanoparticles on a semi-infinite dielectric substrate. In order to better understand the effect of the substrate’s permittivity, ${\varepsilon _2}$, we define the dimensionless parameter as: ${C_i} = \frac{{{\varepsilon _2} - {\varepsilon _1}}}{{{\varepsilon _2} + {\varepsilon _1}}}$. Figure 6(a) shows the variations of the interaction constants versus the contrast factor, ${C_i}$. As mentioned before, the presence of the substrate leads to the creation of non-zero components of interaction constants. Moreover, as the refractive index of the substrate increases, so the magnitude of the image dipole moment approaches the primary dipole. Therefore, the effects of the image dipoles on the interaction constants are increased, as clearly shown in Fig. 6(a).

Next, a parametric study of the effect of the array’s periodicity, on the interaction coefficients, is done for an array of arbitrary nanoparticles located on a glass substrate. It is worth mentioning that, as the image of the reference dipole is not affected by changing the periodicity, we remove its contribution from $\beta ^{\prime}$; and the results are shown in Fig. 6(b). It can be seen from the figure that the interaction constants decrease when the periodicity of the array increases. This is because, by increasing the periodicity of the array, the distance between the adjacent elements increases and they have less electromagnetic influence on each other.

Finally, we investigate the effects of the nanoparticle’s height, on the interaction coefficients. Figure 6(c) shows the result of a 2D periodic array of arbitrary nanoparticles, with a height, *h*, ranging from 100 nm to 200 nm. The array is located on a semi-infinite glass substrate, and the periodicity of the array is 250nm. As shown in Fig. 3, as the nanoparticle’s height increases, so the distance between the primary dipoles and the image dipoles increases. Therefore, the effects of the image dipoles in calculating the interaction fields and thus $\beta ^{\prime}$ are decreased.

After investigation the effect of different parameters on the interaction constants of a substrated metasurface, to analyze the proposed homogenization method for this structure, an array of silicon nano-spheres on a dielectric substrate, with ${\varepsilon _2} = 3{\varepsilon _0}$, is considered. Figure 7 shows the effective polarizabilities extracted by our proposed method in which the individual polarizabilities have been calculated using the method presented in [29]. It is worth highlighting that the analytical Mie theory can only calculate the polarizabilities of special cases such as a spherical nanoparticle in a homogeneous environment; while the far-field method is able to obtain the polarizability tensors of any arbitrary particle — whether in a homogeneous space or on a dielectric substrate. As can be seen from the figure, in contrast to a 2D periodic array of silicon nano-spheres in free space (as see in Fig. 5), which has no electromagnetic coupling, the presence of a substrate breaks the symmetry of the structure and creates an electromagnetic coupling effect. This phenomenon which is referred to as substrate-induced bianisotropy (SIB), has been investigated and reported in different previous works [12,29].

Finally, using Eq. (16) in the Section 2, the transmission and reflection coefficients of the substrated metasurface are calculated. The results are compared with those obtained from a full-wave, Ansys Inc. HFSS simulator; and plotted in Fig. 8(a) which show good agreement. Moreover, we compared the results of the reflection and transmission coefficients with those obtained from free space interaction constants. This comparison is shown in Fig. 8(b). As can be seen from the figure, neglecting the effect of the substrate in the calculation of the interaction constants leads to a larger error when calculating the reflection and transmission coefficients of the substrated metasurface. Therefore, in order to properly analyze the substrated metasurface, it is important to consider the effects of the substrate when calculating the interaction constants.

Characterizing a metasurface by independently considering the effects of individual nanoparticles and the interactions between them, simplify the array analysis. For example, if the array period is modified, only the interaction constants would need to be updated and the other parts remain unchanged. Therefore, the speed of the analysis is higher than other methods as well as commercial software. Moreover, as mentioned before, we used the method presented in [29] to calculate the individual polarizability tensors. Since this method is able to calculate the individual polarizability tensors of any arbitrary dielectric or plasmonic nanoparticles in free space, or on a dielectric substrate, by considering the interaction constants extracted in this paper, one could analyze any arbitrary metasurface located in free space or on a dielectric substrate. Therefore, the proposed approach is general, and these examples are only provided to validate the proposed approach.

In the last example, to examine the capability of the proposed method in dealing with metasurfaces on top of a metallic ground plate, separated by a dielectric spacer, the structure of the metasurface absorber presented in [20], is considered. This structure is an asymmetric Fabry-Perot cavity, with two mirrors; where there is a 2D periodic array of gold patches for the top mirror, and a gold ground plate for the bottom mirror. To analyze this structure, we first, apply the process summarized in the previous section and calculate the reflection and transmission coefficients of the metasurface on an infinite glass substrate. The results are shown in Fig. 9(a). Then, by using the reflection and transmission coefficients calculated with our proposed method, and substituting them into Eq. (17), the reflection coefficient of the structure shown in the inset of Fig. 9(b), can be achieved. The transmission coefficient of this metasurface is completely suppressed due to the sufficiently-thick ground plate. Figure 9(b) also shows the absorption coefficient ($A = 1 - {|R |^2} - {|T |^2}$) of this structure that is calculated using the proposed method, and compared with full wave simulation, HFSS; where a good agreement between the two sets of results is observed.

As discussed in this paper, the proposed method can analyze different substrated metasurfaces using the effect of individual nanoparticles and the interaction between them, independently. Recently, manipulation of the scattering spectrum from metasurfaces has attracted intense attention [18,43–46]. The proposed method in this paper can also be used to synthesis substrated metasurfaces for different applications like perfect absorbers, twist polarizers, asymmetric metasurfaces with different functionalities, and so on, that is planned for our future work.

## 4. Conclusion

In this paper, we have extended the theory of analytical modeling of planar periodic metasurfaces located in free space to a more general case, in the presence of a substrate. In this theory, we have used the advantage of studying the effects of the individual responses of each nanoparticle, as well as the interactions between them, independently. Indeed, we have obtained the new interaction constants that are required for the general analysis of metasurfaces in the presence of a substrate. We focused on substrated metasurfaces with square lattices, and under normal illuminations. Furthermore, as a case study, the effects of the refractive index of substrate, the geometric characteristics of nanoparticles, and the periodicity of the array on the interaction constant tensors, have been analyzed. The validity of the proposed method has been investigated by calculating the characteristic parameters of different metasurfaces; and the reflection and transmission coefficients have been compared to commercial EM solvers, which are in good agreement. We have also demonstrated that our extended theory rather accurately predicts the reflection and transmission coefficients, when compared to the case where we neglect the effects of a substrate in the calculation of the interaction constants. The advantage of the proposed approach is that, by modifying the periodicity of the metasurface, only the interaction constants would need to be modified, and there is no need to repeat the whole analysis process. As a result, the proposed method is faster than commercial software for calculating the transmission and reflection coefficients of metasurfaces with different periodicities. The proposed approach is also general, and can be used for any arbitrary particle in free space or on a dielectric substrate; thus, the examples noted here have only been provided to validate the proposed approach. The analysis of substrated metasurfaces under oblique incidence is planned in our upcoming work.

## Appendix A

Here, we explain the procedure for calculating of other terms in Eq. (6) and also interaction magnetic field created by the image dipoles (${{\bf H}^{int,image}}$).

The re-radiated magnetic field for image of the induced electric dipoles is expressed as follows:

Where ${C_0}$ is the free space light velocity. From the duality theorem, the electric and magnetic fields for image of induced magnetic dipoles are as follows:

By considering ${{\textbf E}_{m^{\prime}}}({\textbf r^{\prime}})$ in Eq. (19) and integration of this field over a sheet with a hole of radius, *R* we have:

Since ${\textbf h}\!=\!{-}h{\textbf z}$ only the parallel components of the magnetic moment in Eq. (20) are considered, which are ${{\textbf m^{\prime}}_\parallel } = {C_i}\begin{array}{{c}} {{{\textbf m}_\parallel }} \end{array}$. Moreover, to calculate ${\textbf E}({{\textbf m^{\prime}}_0},{\textbf h})$ which is the electric field created by the image of the reference magnetic dipole at a distance, *h*, from its center. This field is calculated by substituting ${\textbf n} = {\textbf z}$ and $r^{\prime} = h$ in Eq. (19). The result is:

Similar to the interaction electric field, the interaction magnetic field created by the image dipoles can be explained as follows:

The same procedure can be done for calculating each terms of Eq. (22). The final results are:

As mentioned before, according to the image theorem, ${{\textbf p^{\prime}}_\parallel } ={-} {C_i}\begin{array}{{c}} {{{\textbf p}_\parallel }} \end{array},\begin{array}{{c}} {{{{\textbf p^{\prime}}}_ \bot } = {C_i}\begin{array}{{c}} {{{\textbf p}_ \bot }} \end{array}} \end{array},\begin{array}{{c}} {{{\textbf m}_\parallel } = {C_i}\begin{array}{{c}} {{{\textbf m}_\parallel }} \end{array}} \end{array}$, and ${{\textbf m^{\prime}}_ \bot } ={-} {C_i}\begin{array}{{c}} {{{\textbf m}_ \bot }} \end{array}$.

## Appendix B

To calculate the reflection and transmission coefficients of a metasurface located on a multilayered dielectric substrate, one can extend the multiple reflection method presented in Section 2.2. In this case, as shown in Fig. 10, if additional layers are considered below the dielectric spacer in Fig. 4(a), we need to replace ${r_{23}}$ in Eq. (17) with a generalized reflection coefficient that incorporates sub-surfaces reflections as follows:

Where ${R_{23}}$ is the generalized reflection coefficient that incorporates sub-surfaces reflections. The generalized reflection coefficient for an N layer medium, at the interface between region *i* and region *i+1*, is given by [47]:

Equation (25) is a recursive relation that express ${R_{i, {i + 1} }}$ in terms of ${R_{i + 1, {i + 2} }}$. Moreover, ${r_{i,{i + 1} }}$, ${t_{i, {i + 1} }}$, and ${t_{i + 1, i }}$ are Fresnel reflection coefficients, and they would exist at each boundary if two semi-infinite media from each of the boundary (neglecting the presence of the other boundaries). Using Eq. (25) one can calculate ${R_{23}}$ for a structure with a specified number of layers. Then, by using the reflection and transmission coefficients of a metasurface placed on semi-infinite substrate (${r_{12}}$, ${r_{21}}$, ${t_{12}}$, and ${t_{21}}$) which is described in Section 2.1, the reflection from the metasurface located on a multilayered dielectric substrate are obtained.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **C. Strachan, “The reflexion of light at a surface covered by a monomolecular film,” Math. Proc. Cambridge Philos. Soc. **29**(1), 116–130 (1933). [CrossRef]

**2. **D. V. Sivukhin, “Contribution to the molecular theory of light reflection,” Comptes Rendus (Doklady) Acad. Sci. URSS **36**, 231–234 (1942).

**3. **D. V. Sivukhin, “Molecular theory of the reflection and refraction of light,” Zh. Eksper. Teor. Fiz. **18**, 976–994 (1948).

**4. **L. A. Vainshtein, “On the electrodynamic theory of grids,” Elektronika Bol’shikh Moshchnostei **2**, 14–48 (1963). [CrossRef]

**5. **D. Bedeaux and J. Vlieger, “A statistical theory of the dielectric properties of thin island films: I. The surface material coefficients,” Physica **73**(2), 287–311 (1974). [CrossRef]

**6. **D. C. Langreth, “Macroscopic approach to the theory of surface reflectivity,” Phys. Rev. B **39**(14), 10020–10027 (1989). [CrossRef]

**7. **S. Yoshida, T. Yamaguchi, and A. Kinbara, “Optical properties of aggregated silver films,” J. Opt. Soc. Am. **61**(1), 62–69 (1971). [CrossRef]

**8. **T. Yamaguchi, S. Yoshida, and A. Kinbara, “Optical effect of the substrate on the anomalous absorption of aggregated silver films,” Thin Solid Films **21**(1), 173–187 (1974). [CrossRef]

**9. **C. L. Holloway and E. F. Kuester, “Impedance-type boundary conditions for a periodic interface between a dielectric and a highly conducting medium,” IEEE Trans. Antennas Propag. **48**(10), 1660–1672 (2000). [CrossRef]

**10. **V. A. Kizel, “Modern status of the theory of light reflection,” Sov. Phys. Usp. **10**(4), 485–508 (1968). [CrossRef]

**11. **M. Albooyeh, D. Morits, and C. R. Simovski, “Electromagnetic characterization of substrated metasurfaces,” Metamaterials **5**(4), 178–205 (2011). [CrossRef]

**12. **M. Albooyeh and C. R. Simovski, “Substrate-induced bianisotropy in plasmonic grids,” J. Opt. **13**(10), 105102 (2011). [CrossRef]

**13. **E. F. Kuester, M. A. Mohamed, M. Piket-May, and C. L. Holloway, “Averaged transition conditions for electromagnetic fields at a metafilm,” IEEE Trans. Antennas Propag. **51**(10), 2641–2651 (2003). [CrossRef]

**14. **C. L. Holloway, M. A. Mohamed, E. F. Kuester, and A. Dienstfrey, “Reflection and transmission properties of a metafilm: With an application to a controllable surface composed of resonant particles,” IEEE Trans. Electromagn. Compat. **47**(4), 853–865 (2005). [CrossRef]

**15. **C. L. Holloway, A. Dienstfrey, E. F. Kuester, J. F. O’Hara, A. K. Azad, and A. J. Taylor, “A discussion on the interpretation and characterization of metafilms/metasurfaces: The two-dimensional equivalent of metamaterials,” Metamaterials **3**(2), 100–112 (2009). [CrossRef]

**16. **C. L. Holloway, E. F. Kuester, J. A. Gordon, J. F. O’Hara, J. Booth, and D. R. Smith, “An overview of the theory and applications of metasurfaces: The two-dimensional equivalents of metamaterials,” IEEE Antennas Propag. Mag. **54**(2), 10–35 (2012). [CrossRef]

**17. **M. Yazdi and N. Komjani, “Polarizability calculation of arbitrary individual scatterers, scatterers in arrays, and substrated scatterers,” J. Opt. Soc. Am. **33**(3), 491–500 (2016). [CrossRef]

**18. **Y. Ra’di, V. S. Asadchy, and S. A. Tretyakov, “Total absorption of electromagnetic waves in ultimately thin layers,” IEEE Trans. Antennas Propag. **61**(9), 4606–4614 (2013). [CrossRef]

**19. **M. Yazdi and M. Albooyeh, “Analysis of metasurfaces at oblique incident,” IEEE Trans. Antennas Propag. **65**(5), 2397–2404 (2017). [CrossRef]

**20. **R. Alaee, M. Albooyeh, and C. Rockstuhl, “Theory of metasurface based perfect absorbers,” J. Phys. D: Appl. Phys. **50**(50), 503002 (2017). [CrossRef]

**21. **S. B. Glybovski, S. A. Tretyakov, P. A. Belov, Y. S. Kivshar, and C. R. Simovski, “Metasurface: From microwaves to visible,” Phys. Rep. **634**, 1–72 (2016). [CrossRef]

**22. **N. T. Fofang, T. S. Luk, M. Okandan, G. N. Nielson, and I. Brener, “Substrate-modified scattering properties of silicon nanostructures for solar energy applications,” Opt. Express **21**(4), 4774–4782 (2013). [CrossRef]

**23. **M. Albooyeh and C. R. Simovski, “Huge local field enhancement in perfect plasmonic absorbers,” Opt. Express **20**(20), 21888–21895 (2012). [CrossRef]

**24. **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,” J. Opt. Soc. Am. B **34**(7), D18–D28 (2017). [CrossRef]

**25. **D. Kwon, P. L. Werner, and D. H. Werner, “Optical planar chiral metamaterial designs for strong circular dichroism and polarization rotation,” Opt. Express **16**(16), 11802–11807 (2008). [CrossRef]

**26. **S. Tretyakov, * Analytical modeling in applied electromagnetics* (Artech House, 2003).

**27. **S. I. Maslovski and S. A. Tretyakov, “Full-wave interaction field in two-dimensional arrays of dipole scatterers,” Int. J. Electron. Commun. (AEU) **53**(3), 135–139 (1999).

**28. **C. R. Simovski, M. S. Kondratjev, P. A. Belov, and S. A. Tretyakov, “Interaction effects in two-dimensional bianisotropic arrays,” IEEE Trans. Antennas Propag. **47**(9), 1429–1439 (1999). [CrossRef]

**29. **M. Hesari-Shermeh, B. Abbasi-Arand, and M. Yazdi, “Dipolar analysis of substrated particles using far-field response,” J. Opt. Soc. Am. B **37**(9), 2779–2787 (2020). [CrossRef]

**30. **V. S. Asadchy, I. A. Faniayeu, Y. Radi, and S. A. Tretyakov, “Determining polarizability tensors for an arbitrary small electromagnetic scatterer,” Photonics Nanostruct **12**(4), 298–304 (2014). [CrossRef]

**31. **R. Alaee, M. Albooyeh, M. Yazdi, N. Komjani, C. Simovski, F. Lederer, and C. Rockstuhl, “Magnetoelectric coupling in nonidentical plasmonic nanoparticles: Theory and applications,” Phys. Rev. B **91**(11), 115119 (2015). [CrossRef]

**32. **F. B. Arango and A. F. Koenderink, “Polarizability tensor retrieval for magnetic and plasmonic antenna design,” New J. Phys. **15**(7), 073023 (2013). [CrossRef]

**33. **A. Ishimaru, S. W. Lee, Y. Kuga, and V. Jandhyala, “Generalized constitutive relations for metamaterials based on the Quasi-static Lorentz Theory,” IEEE Trans. Antennas Propag. **51**(10), 2550–2557 (2003). [CrossRef]

**34. **M. S. Mirmoosa, Y. Raadi, V. S. Asadchy, C. R. Simovski, and S. A. Tretyakov, “Polarizabilities of nonreciprocal bianisotropic particles,” Phys. Rev. Appl. **1**(3), 034005 (2014). [CrossRef]

**35. **S. A. Tretyakov, F. Mariotte, C. R. Simovski, T. G. Kharina, and J. P. Heliot, “Analytical antenna model for chiral scatterers: Comparison with numerical and experimental data,” IEEE Trans. Antennas Propag. **44**(7), 1006–1014 (1996). [CrossRef]

**36. **C. R. Simovski, S. A. Tretyakov, A. A. Sochava, B. Sauviac, F. Mariotte, and T. G. Kharina, “Antenna model for conductive omega particles,” J. Electromag. Waves Appl. **11**(11), 1509–1530 (1997). [CrossRef]

**37. **V. Yatsenko and S. Maslovski, “Electromagnetic interaction of parallel arrays of dipole scatterers,” Prog. Electromagn. Res. **25**, 285–307 (2000). [CrossRef]

**38. **P. Yeh, * Optical waves in layered media* (Wiley, 2005).

**39. **H. T. Chen, “Interference theory of metamaterial perfect absorbers,” Opt. Express **20**(7), 7165–7172 (2012). [CrossRef]

**40. **W. T. Doyle, “Optical properties of a suspension of metal spheres,” Phys. Rev. B **39**(14), 9852–9858 (1989). [CrossRef]

**41. **A. B. Evlyukhin, C. Reinhardt, A. Seidel, B. S. Luk’yanchuk, and B. N. Chichkov, “Optical response features of Si-nanoparticle arrays,” Phys. Rev. B **82**(4), 045404 (2010). [CrossRef]

**42. **. Ansys, Inc., “Ansys HFSS 3D Electromagnetic Field Simulator for RF and Wireless Design” (Ansys, Inc., 2020). https://www.ansys.com/products/electronics/ansys-hfss.

**43. **T. Niemi, A. O. Karilainen, and S. A. Tretyakov, “Synthesis of Polarization Transformers,” IEEE Trans. Antennas Propag. **61**(6), 3102–3111 (2013). [CrossRef]

**44. **C. Qu, S. Ma, J. Hao, M. Qiu, X. Li, S. Xiao, Z. Miao, N. Dai, Q. He, S. Sun, and L. Zhou, “Tailor the functionalities of metasurfaces: from perfect absorption to phase modulation,” arXiv preprint arXiv, 1507.00929 (2015).

**45. **H. Guo, J. Lin, M. Qiu, J. Tian, Q. Wang, Y. Li, S. Sun, Q. He, S. Xiao, and L. Zhou, “Flat optical transparent window: mechanism and realization based on metasurfaces,” J. Phys. D: Appl. Phys. **51**(7), 074001 (2018). [CrossRef]

**46. **Y. Li, J. Lin, H. Guo, W. Sun, S. Xiao, and L. Zhou, “A tunable metasurface with switchable functionalities: from perfect transparency to perfect absorption,” Adv. Opt. Mater. **8**(6), 1901548 (2020). [CrossRef]

**47. **W. C. Chew, * Waves and Fields in Inhomogeneous Media* (Van Nostrand Reinhold: New York, 1990).