Nonlinear interaction of two trapped-mode resonances in a bilayer fish-scale metamaterial

We report on a bistable light transmission through a bilayer"fish-scale"(meander-line) metamaterial. It is demonstrated that an all-optical switching may be achieved nearly the frequency of the high-quality-factor Fano-shaped trapped-mode resonance excitation. The nonlinear interaction of two closely spaced trapped-mode resonances in the bilayer structure composed with a Kerr-type nonlinear dielectric slab is analyzed in both frequency and time domains. It is demonstrated that these two resonances react differently on the applied intense light which leads to destination of a multistable transmission.


INTRODUCTION
During recent years research efforts on both fundamental physics and applications of nonlinear optics have demonstrated several major trends [1]. First, the basic concepts of nonlinear optics penetrated into new areas of material science by exploring novel nonlinear materials and nonlinear propagation of light in engineered structures. Second, there are a number of novel concepts related to tunable nonlinear response, and engineered and enhanced nonlinearities, which should be explored more extensively for developing advanced optical tunable nonlinear devices. Moreover, in such structures many of the effects well studied in nonlinear optics can be enhanced by the cavity effects, heterostructures, or Fano-type resonances, that enable much stronger nonlinear response at lower powers. Metamaterials based on planar resonant patterns are expected to become the key building blocks for such nonlinearity-enhanced effects and devices.
It is believed that the fabrication of metal-dielectric composite structures where the dielectric layers possess strong nonlinear response can allow the realization of tunable nonlinear metamaterials. The important area where the role of such nonlinear structures is expected to become more and more pronounced is the creation of all-optical circuitry for computing, information processing, and networking, which is expected to overcome the speed limitations of electronicsbased systems.
Typically metamaterials consist of arrays of magnetically or electrically polarizable elements. In many common configurations, the functional metamaterial inclusions are planar metallic patterns composed of symmetrical split-ring resonators (SRRs), on which currents are induced that flow in response to incident electromagnetic fields. A distinct property of metamaterials is that the capacitive regions of SRRs strongly confine and enhance the local electric field. The capacitive region of a metamaterial element thus provides a natural entry point at which some lumped frequency-dispersive, active, tunable, or nonlinear materials can be introduced, providing a mechanism for coupling the metamaterial geometric resonance to other fundamental material properties [2][3][4]. The hybridization of metamaterials with semiconductors, ferromagnetic/ferrimagnetic materials and other externally tunable materials has already expanded the realm of metamaterials from the passive, linear media initially considered.
While the wave propagation can be controlled by using metamaterials with a unit cell consisting of a single resonant element (symmetrical SRR), other characteristics can also be controlled by using metamaterials composed of coupled resonators. Such metamaterials sometimes are also referred to as electromagnetically induced transparency (EIT)-like metamaterials [5,6]. Typically they consist of identical subwavelength metallic inclusions structured in the form of asymmetrical split-rings (ASRs) [7,8], split squares [9], or their complementary patterns [10]. These elements are arranged periodically and placed on a thin dielectric substrate. It is revealed that since these particles possess specific structural asymmetry, in a certain frequency band, the antiphase current oscillations with almost the same amplitudes appear on the particles parts (arcs). The scattered electromagnetic far field produced by such current oscillations is very weak, which drastically reduces its coupling to free space and therefore radiation losses.
Indeed, both the electric and magnetic dipole radiations of currents oscillating in the arcs of particles are canceled. As a consequence, the strength of the induced current can reach very high value and therefore ensure a high-Q resonant optical response. Such a resonant regime is referred to as so-called "trapped-modes," since this term is traditionally used in describing electromagnetic modes that are weakly coupled to free space.
The features of nonlinear trapped-mode resonances were previously studied theoretically in planar metamaterials composed of ASRs [11] and two concentric rings [12][13][14]. In both configurations the substrate that carries the metallic pattern is considered as a nonlinear dielectric. The effect of nonlinearity appears as the formation of a closed loop of bistable transmission within the frequency of trapped-mode resonance. Since the nonlinear response of these metamaterials operating in the trapped-mode regime is extremely sensitive to the dielectric properties of the substrate it allows one to control switching operation effectively.
Another type of structure, which bears Fano-shaped trapped-mode resonances and is very promising for applications, is planar metamaterials consisting of equidistant array of continuous meander or zigzag metallic strips placed on a thin dielectric substrate (fish-scale structures [15][16][17]). In the past such structures have been investigated both theoretically [15,[18][19][20] and experimentally [15]. It is revealed, when the wave is polarized orthogonally to the strips, the fish-scale structure is transparent across a wide spectral range apart from an isolated wavelength. In contrary, when the wave is polarized parallel to the strips, it becomes a good broadband reflector apart from an isolated wavelength where transmission is high. The same array combined with a homogeneous metallic mirror is also a good reflector which preserves a phase of the reflected wave with respect to the incident wave at this particular wavelength [15,21]. The latter phenomenon is known as a "magnetic mirror" [22]. Finally, the structure also acts as a local field concentrator and a resonant multifold "amplifier" of losses in the constitutive dielectric.
Among others, a trapped-mode resonance can be excited in a fish-scale structure if the incident field is polarized along the strips and when the form of these strips is slightly different from the straight line. In such a system the less the form of strips is different from the straight line, the greater is the quality factor of the trapped-mode resonance. But, unfortunately, as the line curvature decreases, the resonant frequency shifts to the frequency where the Rayleigh anomaly occurs, which diminishes the degree of the field localization within the system.
Nevertheless, in a bilayer configuration of the metamaterial, in which two equal gratings are adjusted to each side of a dielectric slab, it is possible to excite the high-Q trapped-mode resonance far from the Rayleigh anomaly [19]. This resonance appears due to the oscillations of currents that flow in opposite directions relative to the adjacent gratings. In this case the electromagnetic field is strongly localized in the area between the adjacent gratings. Furthermore, in the bilayer structure, besides the trapped-mode resonance excited by gratings, some interference resonances appear in the similar manner as in traditional systems with gratings of straight lines. Therefore, the bilayer fish-scale structure allows one to obtain different resonant features.
The most important thing is that in the bilayer fish-scale structure the trapped-mode resonance has a quality factor that is sufficiently greater than the ones attainable in the convenient interference resonances.
We should also note that the trapped-mode resonances can be excited in a very thin bilayer structure, which is important for practical applications. We argue that the bilayer configuration is of great interest in cases where the substrate is made of a field intensity dependent material (e.g., a Kerr-type medium) because the strong field localization between the gratings can significantly enhance the nonlinear effects. It leads to some particular manifestation of effects of bistability and multistability nearly the frequencies of Fano-shaped trapped-mode resonances, which is presented in this paper.

PROBLEM FORMULATION
A. Bilayer Fish-Scale Metamaterial A studied structure consists of two gratings of planar perfectly conducting infinite wavy-line strips placed on the both sides of a dielectric slab with thickness h and permittivity ε (see Fig. 1). The elementary translation (unit) cell of the structure under study is a square with sides d d x d y . The total length of the strip within the elementary translation cell is S. Suppose that the thickness h and size d are less than the wavelength λ of the incident electromagnetic radiation (h ≪ λ, d < λ). The width of the metal strips and their deviation from the straight line are 2w and Δ, respectively.
Assume that the incident field is a plane monochromatic wave where ⃗p is the unit vector that defines the polarization of the incident wave, jpj 1; A 0 and ⃗ k inc ⃗ e x k inc x ⃗ e y k inc y ⃗ e z k inc z are the magnitude and wave vector of the incident wave, respectively; ⃗ e x , ⃗ e y , ⃗ e z are orts of the coordinate system; k inc k ω ε 0 μ 0 p . The time dependence is supposed to be expiωt.
In the frequency domain a solution of the diffraction problem on the bilayer structure is derived using the method of moments. It is formulated similarly to that one obtained for a single-layer metamaterial's configuration [18,19]. Further, the solution of the linear problem is used to construct a solution of the nonlinear problem.
We shall seek for the field through the entire space in the form of a superposition of the field without gratings ⃗ E d and the field scattered by gratings ⃗ E s , Fig. 1. Fragment of a planar bilayer fish-scale metamaterial and its unit cell. Here, are the electric field strength in the areas above (z ≥ 0) and below (z ≤ −h) the structure, respectively; ⃗ k s fk inc x ; k inc y ; −k inc z g is the wave vector of the reflected wave, and vectors ⃗ R and ⃗ T are defined using corresponding conditions for dielectric boundaries.
Further, we restrict ourself with a case of normal incidence of x-polarized wave (k inc For such a polarization, the reflection and transmission coefficients for a dielectric slab without gratings are The surface current densities on the metal strips in the planes z 0 and z −h are denoted by ⃗ ρ is imposed by the strip midline. In the framework of the method of moments, the metal pattern is treated as a perfect conductor. Based on the boundary conditions of equality to zero of tangential components of the electric field on the strips, we arrive to a system of two integral equations related to the currents induced on the strips: where the linear operators f ⃗ E 1t ⃗ J 1 ; ⃗ J 2 g⃗ κ and f ⃗ E 2t ⃗ J 1 ; ⃗ J 2 g⃗ κ of the linear problem are associated to the Fourier transforms of the surface current density and the tangential components of the radiation field; ⃗ ρ mn md x ⃗ e x nd y ⃗ e y is the radius-vector directed from the unit cell located at the center of coordinates [we shall number it as (0, 0)] to the unit cell located at the position (m, n) for the upper and lower gratings, respectively; ⃗ κ is the wave vector of spatial spectrum of the field.
Using the method of moments, the system of integral Eqs. (3) can be reduced to the system of linear algebraic equations related to the unknown expansion coefficients of the surface currents densities on the strips of the first and second gratings. When these expansion coefficients of the surface currents densities are found (in the same way as it was done for a single-layer metamaterial [18]), it is easy to calculate the magnitude of the spatial harmonics of the field in the form of some plane waves expansion. It can be derived both above and below the structure and If the period of the gratings is less than the wavelength d∕λ < 1, then in the reflected and transmitted fields in the far-field zone there exists only the fundamental spatial harmonics ( ⃗a 00 , ⃗ b 00 ), which propagate along the normal to the structure (single-wave mode). On the other hand, at the chosen polarization of the incident field, the studied structure does not provide any polarization conversion of the scattered field. So, the magnitudes ⃗a 00 and ⃗ b 00 become scalar quantities, and the reflection and transmission coefficients can be defined as R a 00 ∕A 0 and T b 00 ∕A 0 . Consequently, in this paper, the analysis of the diffraction problem is carried out in this single-wave mode (for extra details the reader is referred to [19]).
Thus, the obtained solution allows us to calculate the magnitude and distribution of the current J along the strips, the reflection R and transmission T coefficients as functions of normalized frequency (ae d∕λ), permittivity (ε), and other parameters of the structure: It is noteworthy that the obtained solution has a good experimental confirmation [15].

SPECTRAL PROPERTIES AND BISTABILITY A. Trapped-Mode Resonances
Trapped-mode resonances can be excited in planar metamaterials with particles of different shapes in the case when these particles possess specific structural asymmetry. These resonances are the result of antiphase current oscillations with almost the same magnitudes that appear on the particles parts (arcs). The scattered electromagnetic far field produced by such current oscillations is very weak, which drastically reduces the field coupling to free space and therefore radiation losses. As a consequence, the strength of the induced current can reach very high value and therefore ensure a high-Q resonant optical response. Such a resonant regime is referred to as so-called "trapped-modes." Due to the bilayer configuration of the structure under study there are two possible current distributions that correspond to the trapped-mode resonances. The first distribution is the antiphase current oscillations in arcs of each grating. In this case the structure can be considered as a system of two coupled resonators that operate at the same frequency because the gratings are identical. We have labeled this resonant frequency in Figs. 2 and 3 by the letter ae 1 . Obviously the distance h between the gratings strongly determines the resonant frequency position since this parameter defines the electromagnetic coupling strength. The Q-factor of this resonance is higher in the bilayer structure than the one that existed in a single-layer structure [18] but their similarity lies in the fact that the current magnitude in the metallic pattern depends relatively weakly on the thickness and permittivity of the substrate.
The second distribution is the antiphase current oscillations excited between two adjacent gratings. Similarly we have labeled this resonant frequency in Figs. 2 and 3 by the letter ae 2 . As is known, the closer are the interacting metallic elements, the higher is the Q-factor of the trapped-mode resonance. Thus varying both the distance between the gratings and substrate permittivity changes the trapped-mode resonant conditions and this changing manifests itself in the current magnitude. Remarkably, in the situation of such current distribution, the field is localized between the gratings, i.e., directly in the substrate, which can sufficiently enhance the nonlinear effects if the slab is made of a field intensity dependent material.

B. Nonlinear Spectral Properties
Next we suppose that the structure substrate is a Kerr nonlinear dielectric which permittivity ε depends on the intensity of the electromagnetic field inside it (ε ε 1 ε 2 jE in j 2 ). In this study, we suppose that operations are provided in the infrared part of the spectrum where the trapped-mode resonances in metamaterials with metallic pattern are still well observed and have a high Q-factor [9]. In this range the nonlinear slab can be made of some semiconductor material (e.g., GaAs or InSb). Based on [23], we expect that nonlinear response of the proposed structure becomes apparent at the input light intensity about 1-10 kW cm −2 . Nevertheless, the effects discussed in the article can be actual to a broader range of structures including systems in the microwave range, where nonlinearity is brought about by lumped elements (e.g., transistor-loaded metamaterials).
Under rigorous consideration, the nonlinear permittivity ε of the substrate is inhomogeneous. The permittivity reaches its maximum value directly under the metallic pattern and along the strips this permittivity is also different. Nevertheless at the trapped-mode resonance, the flow of electromagnetic energy is confined to a very small region between the strips, and the crucial impact of the permittivity on the system features occurs in this place. Therefore, it is possible to construct an approximate treatment to estimate the inner-field intensity and to solve the nonlinear problem. This approximate treatment is based on two points. The first one takes into account the transmission line theory and postulates that the inner-field intensity is directly proportional to the square of the current magnitude averaged over a metal pattern extent [12,18]. The second approximation is related to the homogenization procedure which is convenient to the metamaterial methodology [4]. It assumes that, in view of the smallness of the elementary translation cell of the array (d < λ), the nonlinear substrate remains a homogeneous dielectric slab under an action of the intensive light.
In more details, according to this theory, the reference meander wire is considered as a conducting wire periodically loaded by short-circuited sections of transmission lines (each section represents one period). Along these wires the currents flow in opposite directions. Thus, the inner electric field strength is defined as E in V ∕D, where V ZJ is the line voltage, D is the distance between wires in the transmission line,J is the current magnitude averaged along the wire. The impedance is determined at the corresponding resonant frequency ae j d∕λ j , j 1; 2. So, it is Z iZ 0 ∕πdcosh −1 D∕2τ 0 tankΔ, where Z 0 is the impedance of free space, and D d∕2 and D h at the frequencies ae 1 and ae 2 , respectively; τ 0 is the wire radius, and Δ is the length of the equivalent line section. For the equivalent wire radius there is a simple estimate τ 0 w∕2 [24].
From this model it follows that the electric field strength between the wires is directly proportional to the average current magnitudeJ (I in jE in j 2 ∼J 2 ). Thus, the nonlinear equation on the average current magnitude in the metallic pattern can be obtained in the form [11][12][13] J ÃFJ ae; ε 1 ε 2 I in J; whereÃ is a dimensionless coefficient which depicts how many times the incident field magnitude A 0 is greater than 1 V cm −1 . The input field magnitude A 0 is a parameter of this nonlinear equation. So, at a fixed frequency ae, the solution of this equation gives us the averaged current magnitudeJ, which depends on the magnitude of the incident field A 0 . When implementing our computational algorithm, the approach described in Subsection 2.A is included as a part inside the subroutine of the nonlinear Eq. (6). The input parameters of this subroutine are the frequency ae, magnitude A 0 , and inner-field intensity I in . The output parameter is the residual between the inner-field intensity passed to the subroutine and the one calculated inside the subroutine body. The subroutine is called iteratively until the resulting residual becomes less than a predetermined threshold. Note that such a solution may give more than one result, which is a distinctive feature of the nonlinear equation. Thus, on the basis of the currentJA 0 found by a numerical solution of the nonlinear equation, the actual value of permittivity ε of the nonlinear slab is determined and the reflection R and transmission T coefficients are calculated as functions of the frequency ae and the magnitude of the incident field A 0 : One can see that as the magnitude of the incident field rises, the frequency dependences of the inner intensity acquire a form of bent resonances (Fig. 4) and so-called bistable regime occurs. The nature of this effect is studied quite well [11][12][13]25]. Among others, the resonant frequencies of a system are crucially determined by the permittivity of the material of the substrate. Thus, when the frequency of the incident wave is tuned nearly to the resonant frequency, the field localization produces growing the inner light intensity, which can alter the permittivity enough to shift the resonant frequency. When this shift brings the excitation closer to the resonant condition, even more field is localized in the system, which further enhances the shift of the resonance. This positive feedback leads to formation of the hysteresis loop in the inner-field intensity with respect to the incident field magnitude, and, as a result, under a certain magnitude of the incident field, the frequency dependences of the inner-field intensity take a form of bent resonances.
An important point is that this bending is different for distinctive resonances due to difference in their current magnitudes. It results in a specific distortion of the curves of the transmission coefficient magnitude nearly the trapped-mode resonant frequencies. Thus, at the frequency ae 1 ∼ 0.78 the inner field produced by antiphase current oscillations is confined in the area in the vicinity of each grating and it weakly affects the permittivity of the dielectric substrate. In this case the resonant curve acquires a transformation into a closed loop, which is a dedicated characteristic of sharp nonlinear Fano-shaped resonances [26]. The second resonance ae 2 ∼ 0.82 is smooth but the current oscillations produce the strong field concentration between two adjacent gratings directly inside the dielectric substrate. It leads to a considerable distortion of the transmission coefficient magnitude in a wide frequency range, and at a certain incident field magnitude this resonance reaches the first one and tends to overlap it [ Fig. 4(b)]. Evidently, in this case, the transmission coefficient acquires more than two stable states, i.e., the effect of multistability arises.

TIME-DOMAIN SIMULATION AND ALL-OPTICAL SWITCHING
A. Numerical Technique Description The existence of bistability is a necessary but by no means a sufficient condition for practical realization of bistable switching, since the conditions needed to excite a particular branch of the hysteresis loop remain to be determined. Moreover, a hysteresis in the frequency is of limited use because a continuous sweep of the frequency would often be impractical in communication applications. On the other hand, modeling in the frequency domain gives no information on dynamics of response and pulse propagation through a nonlinear metamaterial. Therefore, in this section we demonstrate bistable switching in a bilayer fish-scale metamaterial using timedomain simulations.
On the basis of the microscopic Lorentz-theory approach [27] the effective dielectric permittivity ε eff can be obtained in order to study dynamic behavior of the proposed structure in the time domain (see Appendix A): where ε is permittivity of the nonlinear substrate, and ε 1 is its linear part. Further, we assume the problem parameters in Eq. (8) to be (in the normalized units): ae 1 0.785, ae 2 0.82, α 1 0.01, α 2 0.075, γ 1 0.0008, γ 2 0.005, and ε 1 3. The resulting effective dielectric permittivity of the homogenized layer and corresponding transmission spectra are shown in Fig. 5. One can see that the obtained curves of the transmission spectra are characterized by two asymmetric Fano-shaped resonances, which qualitatively agrees with the results calculated rigorously using the method of moments (see Fig. 3). In order to model light propagation through such an effective medium, we must obtain the polarization corresponding to ε eff . This polarization can be found if we assume that both resonances are in accordance with the oscillations of polarization in the form Fig. 4. Frequency dependences ae d∕λ of (a) the inner field intensity and (b) the transmission coefficient magnitude for different values of the incident field magnitude (nonlinear regime); h∕d 0.2, ε 1 3, ε 2 0.005 cm 2 kW −1 . Other parameters are the same as in Fig. 3.
where E is the electric field strength and ω k are the resonant frequencies. The resulting electric induction is D εE P 1 P 2 . Note that this approach is analogous to that one used in procedures of the Meep package [28]. Assuming P k p k t expiωt and E At expiωt, Eqs. (9) can be rewritten for the amplitudes of polarization p k , where τ ωt, η k 2i γ k , ζ k ω 2 k ε 1 ∕ε − 1 iγ k , and all the values γ k , ω k , and α k are normalized by the light frequency ω. Equation (10) can now be easily discretized, so that the value of polarization at every time instant τ l lΔτ can be calculated by It is worth noting that the values ζ k contain the nonlinear dielectric permittivity ε ε 1 ε 2 jAj 2 and, hence, depend on the strength of the electric field as well.
The propagation of light wave in the effective medium is governed by the wave equation where, as was mentioned above, D contains the contributions of the polarizations P k . Equation (12) can be solved numerically using the finite-difference time-domain technique. We represent the field strength as E At; z expiωt − kz, where ω is a carrier frequency, k ω∕c is the wavenumber, and then solve equations for the complex magnitude At; z. Finally, introducing dimensionless arguments τ ωt and ξ kz, the scheme of calculation of the magnitude values at the mesh points lΔτ; jΔξ is realized as follows: Here the auxiliary values are a 1 1 − iΔτ, a 2 1 iΔτ, g 2Δτ∕Δξ 2 Δτ 2 , and the term responsible for resonant polarizations is The stability of the algorithm is provided by the standard Courant condition between the step intervals of the mesh; in our notation it can be written as Δτ∕Δξ ≤ n, where n is the refractive index. At the edges of the calculation region, we apply the so-called absorbing boundary conditions using total field/scattered field and perfectly matched layer techniques [29].

B. Dynamics of Response and Pulse Propagation
The results of calculations using the method described above are shown in Fig. 6. The parameters of the medium are the same as previously, except the thickness, which is now taken to be h∕d 0.5. The dimensionless frequency of the incident light is ae 0.783, i.e., it is chosen to be nearly the frequency of the first trapped-mode resonance (Fig. 3). The normalized magnitude A∕A 0 of light varies from 1 to 11 with a certain period which is sufficient for steady-state establishment. Recall that the main feature of the considered fish-scale metamaterial is the appearance of two particular trappedmode resonances. The same two resonances can be seen if we plot the effective permittivity ε eff as a function of the nonlinear contribution ε and, hence, intensity. Therefore, if the range of magnitude variation is large enough, one can expect to obtain two hysteresis loops, which is the indication of multistability. In order to reduce absorption and the necessary input intensities, we took the thinner layer than in the previous example, h∕d 0.2, and get closer to the resonance ae 0.784; other parameters are left the same. The results of the calculations are shown in Fig. 7. It turned out that the appearance of the second loop strongly depends on the initial light magnitude. Starting from low input intensities and increasing it with a certain step, we obtained the transition between distinct stable states corresponding to the shift of the first resonance, but we were not able to observe the similar jump associated with the second resonance [see Fig. 7(a)]. Perhaps the transition magnitude is very high, so that we could not reach it. However, we can accept the opposite strategy: start just from the high enough magnitude and decrease it. In this case, we obtain clear evidence of the second loop, as can be seen in Fig. 7(b). The characteristics of the first loop are identical in both cases.
Another way to study bistable response of our nonlinear metamaterial is to launch a pulse into it and analyze the shape of the resulting radiation profile. We consider a Gaussian pulse with magnitude A A p exp−t 2 ∕2t 2 p , where t p Nλ∕c is its duration measured through the number of periods N. In our calculations, N 200 and A p 20A 0 . The dynamics of transmitted and reflected light are presented in Fig. 8(a). The characteristic features of pulse interaction with bistable medium are seen-the flattening of the profile and the overshoot at the leading edge of the pulse [25]. The bistable switching is not evident in these figures and appears through the asymmetric form of the output pulse. To reconstruct the hysteresis loop, we calculate the ratio of the transmitted intensity to the input one at every time instant (the propagation time through the layer is neglected), and obtain the dependence of the transmission on the input intensity for both rising and decaying pulse slopes. The bistable curve derived from the transmitted pulse profile [ Fig. 8(b)] is in agreement with the results of stationary calculations of Fig. 7(a). To prove the  existence of the second loop, we launched not a full pulse, but only a half of it. The profile of transmitted radiation given in Fig. 8(c) clearly shows the occurrence of two different stable states. The corresponding loop is presented in Fig. 8(d) where the curve for ascending intensities is taken from the panel (b). Thus, the results of stationary and pulsed calculations are in full accordance with each other. Finally, we should touch on the problem of modulation instability (MI) [30][31][32]. Though there is not any evidence of self-pulsations or other MI effects in the stationary regime (Figs. 6 and 7) and in the case of the full pulse [ Fig. 8(a)], the beatings in Fig. 8(c) can be considered as a probable manifestation of MI. To make the final decision, the careful study of the conditions for MI appearance in the structure considered should be done, which we plan to perform in our future investigations.

CONCLUSIONS
The light transmission through a bilayer fish-scale metamaterial is considered in both frequency and time domains. An important feature of the studied structure is its ability to support two distinct configurations of the trapped-mode resonance excitation. The first configuration corresponds to a specific distribution of currents, which oscillate in the same manner on each grating. The second configuration is a result of antiphase oscillations of currents, which flow in opposite directions relative to the adjacent gratings. In the latter case it is found out that the electromagnetic field is confined between the gratings, i.e., directly in the substrate which leads to manifestation of the effects of bistability and multistability. They appear in the form of a significant distortion of transmission coefficient magnitude over a wide frequency band. In addition, a study of the pulse propagation through the metamaterial is investigated. The results of stationary and pulsed calculations are in full compliance with each other.
We argue that more parameters can be controlled in the studied metamaterial and that it has better performance than metamaterials with a unit cell of a single symmetrical resonant element. This independent control provides an alternative route to the optimization of nonlinear materials, which form the basis of optical devices such as mixers, frequency doublers, and modulators.

APPENDIX A
An electron in the metal is affected by the external electric field E, as well as by the Coulomb force F C from the other electric charges q j and electromotive force F L from the electric currents. The Coulomb force is the restoring force described as in the case of harmonic oscillator by F C P j eq j ∕εr j −mω 2 0 ξ, where e and m are the electron charge and mass, ξ is the electron's displacement. It is important that the force and the oscillator (resonant) frequency ω 2 0 are inversely proportional to the dielectric permittivity of the nonlinear substrate ε. In the linear and nonlinear regimes, we have ω 2 01 a∕ε 1 and ω 2 0 a∕ε, respectively, where ε 1 is the linear part of ε. Thus, ω 2 0 ω 2 01 ε 1 ∕ε. The electromotive force is determined by the magnetic field induced by the electric current j eN∂ξ∕∂t and takes the form F L −e∂Φ∕c∂t ∼ ∂B∕∂t ∼ ∂j∕∂t ∼ ∂ 2 ξ∕∂t 2 , where N is the density of electrons, Φ is the magnetic flux. We present this force by means of the coefficient β, so that F L −βm∂ 2 ξ∕∂t 2 −βmξ. Thus, the equation of motion for the electron is mξ mγ _ ξ eE − mω 2 0 ξ − βmξ: (A1) Implying the stationary time dependence expiωt for the field and displacement, we get