## Abstract

A simple semi-analytical single mode model describing the mechanism of metal-insulator-metal gap plasmon mode excitation by a plane wave is proposed. The role of the insulator-metal-insulator (IMI) lattice mode is highlighted. Although many other works addressed this issue, the crucial role of this mode has never been demonstrated before. This mode appears as the missing link that ensures energy transfer between the incident plane wave and the metal-insulator-metal (MIM) gap plasmon mode. In this single mode model, the grating layer, the host layer of the IMI lattice mode, is viewed as an effective homogeneous medium, and the scattering matrix characterizing the interaction between the IMI lattice mode and the MIM gap plasmon mode is easily computed. The proposed simplified model allows us to grasp the physical origin of the modes of the system and to predict accurately the resonance frequencies of the 1D structure. These modes are classified in symmetric and antisymmetric modes. The incident field, by its symmetry properties, acts on the system as a selection rule, activating a class of modes with the same symmetry properties as itself.

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

## 1. Introduction

Metasurfaces consist of periodical subwavelength-sized arrangements of plasmonic resonators. The interaction between an incoming electromagnetic wave and such surfaces, allows for a powerful and precise control of the wavefront and leads to some artificial properties of these surfaces [1]. The electromagnetic light absorption [2] is one of these properties that plays a crucial role in photonics technology. There is currently a great interest in developing artificial and ultra-thin optical material allowing perfect absorption of light on a very large angular aperture [3]. In visible and infrared ranges, plasmonic absorbers have proven to be effective and often involve resonance phenomena in very small gaps, like in the case of gap plasmon resonators [4, 5] or bow-tie antennas [6]. The nanopatch metasurface is one of the most studied and simplest structures that exhibit unique absorption and field enhancement. This structure consists of planar nano-resonators separated from metallic substrate film by an insulator [5,7–10]. Although several studies have been devoted to an explanation and an interpretation of the coupling mechanism between the incident plane wave and the gap plasmon, some aspects of this mechanism are still unclear. Recently, based on a Fabry-Perot cavity model, C. Lemaitre proposed a semi-analytical model to explain the role of the angle of incidence in the excitation of the resonance frequencies of the system. Through their cavity model, C. Lemaitre et al. [11] suggested that two kinds of resonances of the gap plasmon cavity can be distinguished: the odd resonances with a very high absorption cross-section which decreases when the incidence angle increases, and the even resonances that cannot be excited at normal incidence. A year later, X. Jia et al. [12] provided another analytical model to explain the angle-dependence of the absorption resonance of a second, sharp absorption dip of a nanopatch metasurface. They suggested that this second resonance is a collective effect involving the excitation of surface plasmon modes and relates to a Wood’s anomaly. From my point of view, both works [11], [12] are complementary since the second angle-dependent absorption discussed in [12] is nothing more than the excitation of a high-order cavity mode of the Fabry-Perot cavity highlighted in [11]. Put side by side both works constitute a major advance in the attempt to explain the coupling between the plane wave and the gap plasmon mode. However, it is possible to deepen and complete them by a modal analysis of the structure. Thus I present in this paper another interpretation of this phenomenon which relies on the existence and the excitation of a super mode living in the corrugated zone. This insulator-metal-insulator (IMI) lattice mode ensures the energy transfer between the incident field and the metal-insulator-metal (MIM) gap plasmon mode. Since the feature of devices at hand is linked to a plasmon resonance phenomenon, it may be successfully treated as an eigenvalue problem with specific boundary conditions i.e. a boundary value problem. This boundary value problem is efficiently solved throughout a full-wave polynomial modal method (PMM) [13–16]. First, I show, through the modal analysis based on the PMM, that the reflection spectrum of the IMI lattice mode matches very well with the reflectance of the structure. Second, a new physical analysis and interpretation of the reflection spectrum of the system is presented. The proposed model is based on the scattering matrix analysis obtained from the coupling between the specific IMI mode living in an equivalent homogeneous medium and the MIM gap plasmon mode. The existence of two MIM cavity gap plasmon modes classes under the ribbon is outlined: symmetrical and anti-symmetrical cavity gap plasmons and each class of modes satisfies to a specific dispersion relation. I show that the existence of these two classes of modes is intrinsic to the system and therefore does not strongly depend on the angle of incidence. However, the excitation or not of the anti-symmetrical mode is strongly related to the symmetry properties of the incident field.

## 2. Modal analysis of the IMI-MIM coupling

The structure under consideration is presented in figure 1(a). It consists of a subwavelength periodically distributed *a*-width silver nanoribbons with relative permittivity *ε*^{(m1)}, deposited on a *ε*^{(2)}-relative permittivity dielectric spacer layer, and a *ε*^{(m2)} relative permittivity thick gold-ground substrate. The height *h* 1 of the silver-nanoribbons is set to 75*nm* and the gold-ground plane height *h*_{3} = 50*nm*. The spacer layer height *h*_{2} is equal to 5*nm* and its relative permittivity *ε*^{(2)} is set to 1.54^{2}. The relative permittivity of the dielectric material of the grating is denoted by *ε*^{(d)}. In this manuscript, the dispersive relative permittivity functions *ε*^{(m1)} and *ε*^{(m2)} of the metals are described by the Drude-Lorentz model. See references [17–19] for more details. This structure is shined, from the upper medium (with relative permittivity *ε*^{(0)}) by a TM polarized plane wave (**H** = *H _{y}*

**e**). The wave vector of the incident wave is denoted by

_{y}**K**=

_{0}*k*

_{0}(

*α*

_{0}

**e**+

_{x}*β*

_{0}

**e**+

_{y}*γ*

_{0}

**e**) where

_{z}*k*

_{0}= 2

*π*/

*λ*=

*ω*/

*c*is the wavenumber,

*λ*being the wavelength and

*c*the light velocity in vacuum. The reflection of the structure is computed throughout the polynomial modal method (PMM). The structure is divided into 3 layers ${I}_{z}^{(k)}$,

*k*= 1 : 3 in (

*O*,

*z*) direction. In each layer ${I}_{z}^{(k)}$, encapsulated between the planes

*z*=

*z*

_{k−1}and

*z*=

*z*, the general solution

_{k}*ϕ*

^{(k)}(

*x*,

*z*) representing the

*H*component

_{y}- backward waves propagating along decreasing values of
*z*$${H}_{y}^{(k)-}(x,z)=\sum _{q}{B}_{q}^{(k)}{e}^{i{k}_{0}{\gamma}_{q}^{(k)}(z-{z}_{k-1})}\sum _{n}{H}_{nq}^{(k)}{P}_{n}(x)=\sum _{q}{B}_{q}^{(k)}{e}^{i{k}_{0}{\gamma}_{q}^{(k)}(z-{z}_{k-1})}{H}_{q}^{(k)}(x),$$where*P*is a set of polynomial functions allowing to expand the eigenfunctions ${H}_{q}^{(k)}$. These eigenfunctions are computed as eigenvectors of the TM propagation equation: where_{n}

We solve the eigenvalue equation (4) in each layer ${I}_{z}^{(k)}$, *k* = 1 : 3 before writing the boundary conditions in the (*Oz*) direction. These boundary conditions equations are linked thanks to the scattering matrix algorithm. The weighting coefficients ${A}_{q}^{(k)}$ and ${B}_{q}^{(k)}$ in each layer can be computed, and we are able to know which mode is excited in each layer. At this stage, we assert that in the grating zone, a particular mode mainly ensures the energy transfer between the incident plane wave and the MIM gap plasmon mode under the ribbons. As a starting point of our analysis, we consider the most slowly decaying evanescent mode of layer ${I}_{z}^{(1)}$ whose effective index is denoted by ${\gamma}_{0}^{(1)}$. As shows Fig. (2(b)), this mode decreases in the silver metal and strongly depends on the parameter *d*/*λ* in the space between the ribbons. We define in Eq. (6) below, at the bottom interface of layer ${I}_{z}^{(1)}$, *i.e.* at *z* = *z*_{1}, a coefficient of reflection ${r}_{0}^{(1)}$ (see Fig. (1(a))) of the eigenmode which effective index is ${\gamma}_{0}^{(1)}$:

*λ*∈ [0.6, 1.5]

*μm*). Comparing ${\left|{r}_{0}^{(1)}\right|}^{2}$ with the reflectance of the structure, we remark that the spectral response of this lattice IMI mode ${\gamma}_{0}^{(1)}$ accurately matches with the reflectance of the whole system. The IMI mode is responsible of the Fabry-Perot cavity input and output ports excitation. In the following section we give a simple interpretation of the modes of the structure in terms of an approximated coupling scheme and we compare our model with the exact numerical results obtained from the PMM.

## 3. Study of the coupling between the IMI lattice mode and the MIM gap plasmon mode

In the current case, the transversal geometrical parameters of the grating are smaller than the incident field wavelength *λ* (*d* << *λ*). Therefore the host medium of ${\gamma}_{0}^{(1)}$ IMI lattice mode can be approximated by a homogeneous medium with equivalent permittivity *ε*^{(1)} = 〈1/*ε*^{(m1,d)}(*x*)〉^{−1} [19]. It is then possible to define for this mode an effective index ${\alpha}_{0}^{(1)}$ along (*Ox*)-axis as follows:

*a*-width

*ε*

^{(2)}-permittivity media where lives the ${\alpha}_{0}^{(2)}$-effective index gap plasmon mode [16], is encapsulated between an input and output

*ε*

^{(1)}-permittivity media. The scattering matrix parameters of this system, defined in terms of incident and reflected waves at input and output ports, are written as follows

*ϕ*=

*ϕ*. The phase ${\varphi}_{a}={e}^{-i{k}_{0}{\alpha}_{0}^{(2)}a}$ handles the propagation of the MIM gap plasmon mode with effective index ${\alpha}_{0}^{(2)}$ under the

_{a}ϕ_{b}*a*-width patch line. As it is shown in Fig. (1(b)), the phase references of the coefficients

*C*and

_{k}*D*,

_{k}*k*= 1, 2, are shifted away to the limits of the

*ε*

^{(2)}-slab. Therefore we add an additional phase delay

*ϕ*at each interface of this slab. For a given angle of incidence, this phase delay

_{b}*ϕ*, that will be computed later, mainly depends on the transversal effective index ${\alpha}_{0}^{(1)}$ of the IMI lattice mode, the width

_{b}*b*between two ribbons and on the period

*d*. The coefficients

*r*and

_{i}*t*,

_{i}*i*= 1, 2 are the Fresnel coefficients at the interface

*ε*

^{(i)}/

*ε*

^{(i+1)}in TM polarization, ${r}_{i}=\frac{{\alpha}_{0}^{(i)}/{\epsilon}^{(i)}-{\alpha}_{0}^{(i+1)}/{\epsilon}^{(i+1)}}{{\alpha}_{0}^{(i)}/{\epsilon}^{(i)}+{\alpha}_{0}^{(i+1)}/{\epsilon}^{(i+1)}}$,${t}_{i}=\frac{2{\alpha}_{0}^{(i)}/{\epsilon}^{(i)}}{{\alpha}_{0}^{(i)}/{\epsilon}^{(i)}+{\alpha}_{0}^{(i+1)}/{\epsilon}^{(i+1)}}$. In this paper,

*ε*

^{(3)}=

*ε*

^{(1)}and ${\alpha}_{0}^{(3)}={\alpha}_{0}^{(1)}$. Now, let us set

*S*-matrix parameters are then written in the following simpler manner

*r*

_{1}(

*ω*) = (1 −

*n*(

*ω*))/(1 +

*n*(

*ω*)),

*t*

_{1}= 2/(1 +

*n*(

*ω*)),

*t*

_{2}=

*n*(

*ω*)

*t*

_{1};

*ω*) of the matrix-valued function

*S*(

*ω*) of equation Eq. (8):

*ϕ*is required. As said before, for a given incident wave, it is easy to conceive that this phase depends on the parameters

_{b}*b*as well as on the period

*d*. We suggest, the following analytical form: where

*f*(

*θ*

_{0},

*d*) is estimated through a linear interpolation method, that generally, takes two data points, say (

*d*

_{1},

*f*(

*θ*

_{0},

*d*

_{1}) and (

*d*

_{2},

*f*(

*θ*

_{0},

*d*

_{2}). The interpolant is then given by:

*f*(

*θ*

_{0}= 0

^{0},

*d*

_{1}= 200

*nm*) = 1 and

*f*(

*θ*

_{0}= 0

^{0},

*d*

_{2}= 300

*nm*) =

*b*/

*d*

_{2}. Figures (3(a)) and (3(b)) show the spectra of |Δ(

*λ*)| = |

*S*

_{11}(

*λ*)

*S*

_{22}(

*λ*)−

*S*

_{12}(

*λ*)

*S*

_{21}(

*λ*)| (with and without the additional phase

*ϕ*), |

_{b}*S*

_{11}(

*λ*)−

*S*

_{12}(

*λ*)| and |

*S*

_{11}(

*λ*) +

*S*

_{12}(

*λ*)|. The reflectance of the structure computed with the PMM is also reported in these figures. For these numerical simulations, we consider a period

*d*of 0.25

*μm*, the incident angle is set to

*θ*

_{0}= 0° in Fig. (3(a)) and

*θ*

_{0}= 40° in Fig. (3(b)).

First, we remark that:

- whatever the incidence, two resonance wavelengths are predicted by Eq. (12): a first one close to
*λ*= 1.01*μm*and a second one close to*λ*= 0.672*μm* - each resonance wavelength of the structure coincides with a minimum of the function |
*S*_{11}−*S*_{12}| or |*S*_{11}+*S*_{12}| - however at the normal incidence (
*θ*_{0}= 0^{0}) the reflectance of the structure exhibits only one resonance phenomenon close to*λ*= 1.01*μm*while for*θ*_{0}= 40^{0}, both predicted resonances phenomena coincide with the minimums of the reflection curve.

Before explaining each point raised above, let us first begin by clarifying and classifying the different solutions of Eqs.(12). Equations (12) exhibit two kinds of solution denoted *ϕ*^{+} and *ϕ*^{−} satisfying:

*x*. The solution

*ϕ*

^{+}(

*ω*)

*r*

_{1}(

*ω*) = ±1 corresponds to a

*x*-direction forward wave while the solution [

*ϕ*

^{−}(

*ω*)]

^{−1}

*r*

_{1}(

*ω*) = ±1 corresponds to a backward wave solution. Among these both solutions, we will focus on the forward wave

*ϕ*

^{+}(

*ω*) since both solutions are equivalent. Considering Eqs (12) and (15), one can distinguish two classes of solution

*ϕ*

^{+}(

*ω*):

As it is shown in Figs. (3(a)) and (3(b)), each resonance frequency corresponds to a solution *ϕ*^{s+} or *ϕ*^{a+}. The solutions *ϕ*^{s+} of *S*_{11}(*ω*) − *S*_{12}(*ω*) = 0 correspond to a configuration in which the output and input ports of the slab are excited by two fields of amplitudes *C*_{1} and *C*_{2} with a phase shift of (2*p* + 1)*π*, *p* ∈ . In this case, the field in the cavity has a symmetrical shape with respect to the *z*-median plane of the nanoribbon. Similarly, the relation *S*_{11}(*ω*) + *S*_{12}(*ω*) = 0 leads to a phase delay of 2*pπ*, (*p* ∈ ) between both ports of the slab. In this second case, there is in the MIM gap, an anti-symmetrical cavity mode, with respect to the *z*-median plane of the ribbon. Therefore

*ϕ*

^{s+}corresponds to a symmetrical gap plasmon cavity mode while

*ϕ*

^{a+}is an anti-symmetrical gap plasmon cavity mode. We support our assertion with Figs (3(c)) and (3(d)) where are plotted, the real part of magnetic field

*H*(

_{y}*x*,

*z*) corresponding to resonance wavelengths

*λ*= 1.01

*μm*and

*λ*= 0.67

*μm*. The incidence angle is set to 40°. For

*λ*= 1.01

*μm*a symmetrical gap plasmon cavity mode is excited by the incident plane wave while an anti-symmetrical gap plasmon cavity mode is excited by the incident plane wave at

*λ*= 0.67

*μm*. It thus becomes obvious that the resonance wavelength close to

*λ*= 0.67

*μm*corresponding to an anti-symmetrical solutions

*ϕ*

^{a+}cannot be excited from a symmetrical incident wave such as an incident plane wave at normal incidence.

Note that the additional phase *ϕ*_{b} does not have a major impact on the existence of the solutions of equation Eq. (11). As shown in Fig. (3(a)) and (3(b)), if *ϕ _{b}* is omitted in the approximate model

*i.e. ϕ*is set to 1, the positions of the resonance frequencies are only shifted. The existence of the solutions

_{b}*ϕ*

^{s+}and

*ϕ*

^{a+}is not related to the angle of incidence. The incidence angle only and mainly modifies the positions of these resonance frequencies. To summarize, equation (11) provides a necessary condition for the Δ(

*ω*) function being minimum, and hence for all possible resonance frequencies of the system. This condition is based on the knowledge of the

*S*-matrix resulting from the interaction between the IMI lattice mode and MIM gap plasmon mode alone. The condition of Eq. (11) predicts the resonance frequencies but does not handle the energy quantity conveyed at these frequencies. The incident field acts on the system as a selection rule. Secondly, considering again Figs. (3(a)) and (3(b)). One can also remark that the function doesn’t really vanish at its minimum. We assert that the function |Δ(

*ω*)| doesn’t really vanish at its minimum because the current model does not take into account supplementary losses due to the

*ε*

^{(m2)}metal,

*i.e.*gold. To prove this assertion, we introduce artificial losses through the wavenumber

*k*0 appearing in the phases

*ϕ*and

_{a}*ϕ*: (4(a))

_{b}*λ*,

*k″*

_{0}) for

*λ*∈ [0.6, 3]

*μm*and

*k″*

_{0}∈ [0, 1]. Regarding this figure, and for resonance wavelength

*λ*= 1.01

*μm*, the modulus of Δ decreases to a minimum value very close to zero. As reported in Fig. (4(b)), the optimal value of

*k″*

_{0}for the current example is

*k″*

_{0}= 0.2

*rad*/

*μm*.

Finally, by using our model, we provide some numerical simulations in a different cases by modifying the period *d*, and the parameter *a*. See Figs. (5(a)) and (5(b)). In these figures, we compare the spectra of |Δ(*λ*)|, |*S*_{11}(*λ*) − *S*_{12}(*λ*)|, |*S*_{11}(*λ*) + *S*_{12}(*λ*)| with the reflection curve obtained from PMM (solid line), for the following numerical parameters: *d* = 300*nm*, *a* = 150*nm*. The incidence angle is set to 0°, in Fig. (5(a)) and 40° in Fig. (5(b)). For *λ* ∈ [0.6, 3], three resonance wavelengths can be distinguished: *λ*_{1} ≃ 1.736*μm* (point 1), *λ*_{2} ≃ 0.97*μm* (point 2) and *λ*_{3} ≃ 0.752*μm* (point 3). All these resonance wavelengths are accurately predicted by our single mode model. As it is shown in Figs. (5(c)) and (5(d)), the point 2 corresponds to an anti-symmetrical cavity mode resonance while at point 3 a symmetrical cavity gap-plasmon mode is excited. As it is expected, only the symmetrical modes are excited at normal incidence. See Fig. (5(a)).

## 4. Conclusion

In conclusion, we have proposed a simple model, providing another interpretation of the mechanism of metal-insulator-metal gap plasmon mode resonance of a nanoribbons periodic array. In the proposed model, the role of insulator-metal-insulator lattice mode is highlighted. We shown through a rigorous numerical modal method that this mode appears as the main channel allowing the energy transfer from the plane wave to the gap plasmon mode under the ribbons. We then proposed a single mode model allowing to confirm the role of the IMI mode. We have shown that the resonance frequencies of the structure can be classified into two classes, namely: symmetrical and anti-symmetrical modes which existence may be predicted independently from the incident field parameters. The incident field acts as a selection rule allowing exaltation or annihilation of such modes. A symmetrical incident field (field at normal incidence for example), cannot excite anti-symmetrical modes. The proposed simple model is reasonably a good approximation that allows to grasp the underlying physical phenomena occurring in the system. This model could be extended to 2D periodic arrays of resonators with arbitrary shape through a modal method implemented in the matched coordinates system [20].

## Funding

Agence Nationale de la Recherche.

## Acknowledgments

This work has been sponsored by the French government research program “Investissements d’Avenir” through the IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25)

## References

**1. **N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nat. Mater. **13**, 139–150 (2014). [CrossRef] [PubMed]

**2. **A. Hartstein, J.R. Kirtley, and J.C. Tsang, “Enhancement in the infrared absorption from molecular monolayers with thin metal overlayers,” Phys. Rev. Lett. **45**, 201–204 (1980) [CrossRef]

**3. **R. Smaali, F. Omeis, A. Moreau, T. Taliercio, and E. Centeno, “A universal design to realize a tunable perfect absorber from infrared to microwaves, ” Sci. Rep. **6**, 32589 (2016). [CrossRef] [PubMed]

**4. **M. G. Nielsen, D. K. Gramotnev, A. Pors, O. Albreksten, and S. I. Bozhevolnyi, “Continuous layer gap plasmon resonators,” Opt. Express **19**, 19310–19322 (2011). [CrossRef] [PubMed]

**5. **A. Moreau, C. Ciracì, J. J. Mock, R. T. Hill, Q. Wang, B. J. Wiley, A. Chilkoti, and D. R. Smith, “Controlled-reflectance surfaces with film-coupled colloidal nanoantennas,” Nature **492**, 86–89 (2012). [CrossRef] [PubMed]

**6. **N. Yu, E. Cubukcu, L. Diehl, D. Bour, S. Corzine, J. Zhu, G. Hőfler, K. B. Crozier, and F. Capasso, “Bowtie plasmonic quantum cascade laser antenna,” Opt. Express **15**, 13272–13281 (2007). [CrossRef] [PubMed]

**7. **T. Sondergaard and S. Bozhevolnyi, “Slow-plasmon resonant nanostructures: Scattering and field enhancements,” Phys. Rev. B **75**, 073402 (2007). [CrossRef]

**8. **J. Yang, C. Sauvan, A. Jouanin, S. Collin, J.-L. Pelouard, and P. Lalanne, “Ultra small metal-insulator-metal nano resonators: impact of slow-wave effects on the quality factor,” Opt. Express **20**, 16880–16891 (2012). [CrossRef]

**9. **R. Esteban, T. V. Teperik, and J. J. Greffet, “Optical Patch Antennas for Single Photon Emission Using Surface Plasmon Resonances,” Phys. Rev. Lett. **104**, 026802 (2010). [CrossRef] [PubMed]

**10. **C. Ciraci, J. B. Lassiter, A. Moreau, and D. R. Smith, “Quasi-analytic study of scattering from optical plasmonic patch antennas,” J. Appl. Phys. **114**, 163108 (2013). [CrossRef]

**11. **C. Lemaitre, E. Centeno, and A. Moreau, “Interferometric control of the absorption in optical patch antennas,” Sci. Rep. **26**, 3004–3012 (2017).

**12. **X. Jia, P. Bowen, Z. Huang, X. Liu, C. Bingham, and D. R. Smith, “Clarification of surface modes of a periodic nanopatch metasurface,” Opt. Express **26**, 3004–3012 (2018). [CrossRef] [PubMed]

**13. **K. Edee, “Modal method based on subsectional Gegenbauer polynomial expansion for lamellar gratings,” J. Opt. Soc. Am. A **28**, 2006–2013 (2011). [CrossRef]

**14. **K. Edee, I. Fenniche, G. Granet, and B. Guizal, “Modal method based on subsectional Gegenbauer polynomial expansion for lamellar gratings: weighting function, convergence and stability,” Prog. Electromagn. Res. **133**, 17–35 (2012). [CrossRef]

**15. **K. Edee and J.-P. Plumey, “Numerical scheme for the modal method based on subsectional Gegenbauer polynomial expansion: application to biperiodic binary grating,” J. Opt. Soc. Am. A **32**, 402–410 (2015). [CrossRef]

**16. **K. Edee, “Polynomial modal method for analysis of the coupling between a gap plasmon waveguide and a square ring resonator,” J. Appl. Phys. **122**, 153102 (2017). [CrossRef]

**17. **R. Brendel and D. Bormann, “An infrared dielectric function model for amorphous solids,” J. Appl. Phys **71**, 1–6 (1992). [CrossRef]

**18. **A. D. Rakiá, A. B. Djurišić, J. M. Elazar, and M. L. Majewski, “Optical properties of metallic films for vertical-cavity optoelectronic devices,” Appl. Opt. **37**, 5271–5283 (1998). [CrossRef]

**19. **K. Edee, “Single mode approach with versatile surface wave phase correction for the extraordinary optical transmission comprehension of 1D period nano-slits arrays,” Opt. Soc. Am. Cont. **1**, 613–624 (2018).

**20. **K. Edee, J.-P. Plumey, A. Moreau, and B. Guizal, “Matched coordinates in the framework of polynomial modal methods for complex metasurfaces modeling,” J. Opt. Soc. Am. A **35**, 608–615 (2018). [CrossRef]