{'fraction_non_alphanumeric': 0.027783591460862286, 'fraction_numerical': 0.011458769359564672, 'mean_word_length': 4.527183342972816, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dispersion Properties of Multilayered Metal-Dieletric Metamaterials

We study dispersion properties of layered metal-dielectric media having different layers thicknesses ratios. Plotting dispersion diagrams and isofrequency contours, we find that strong nonlocality is inherent property of such periodic structures. We introduce different level of losses and analyze complex modes of the metamaterial demonstrating that it operates in a regime with infinite numbers of eigenmodes, with nonlocality slightly affected by losses.

It is known that the structure under consideration demonstrates a strong nonlocal response [19]. Previously, we have shown in [20] that under certain conditions two extraordinary waves instead of one can be excited in the metamaterial by incident TM-polarized light beam. The existence of such an additional light wave is a clear manifestation of the presence of spatial dispersion in the metamaterial, and it validates strong nonlocal nature of the structure. In the recent paper [21] this phenomenon was treated in the framework of arising field of nonlocal transformation optics. An important question of will the nonlocality and the beam splitting phenomenon still remain in real structures with introduction of losses emerges though. Consequently, it is needed to study complex modes of the multilayered metamaterial and clarify an effect of losses. Local effective medium model along with transfer matrix method [22] are used in order to describe the structure. The former models the structure depicted in Fig. 1 with layers permittivities ε 1 , ε 2 and thicknesses d 1 , d 2 as the uniaxial anisotropic crystal of the following form: In the case when ε and ε ⊥ are of opposite sign, i.e. ε ε ⊥ < 0, a material is called indefinite [1] or hyperbolic [23,24]. Since hyperbolic metamaterials allow high-k modes they may have very high photonic density of states [25,26].
However, it should be noticed that the local effective medium model is not accurate and used in the present paper only to show nonlocal behaviour of the multilayered metamaterial by comparison between results provided by the local model and the actual description. The nonlocal effective medium model (see Ref. [27]) is capable of providing analytic results coinciding with ones of transfer matrix method which is the most correct electromagnetic description of multilayers.
Analysis of MMDN was accomplished in the work [28] for lossless metamaterials. However, proper attention was not devoted to the imprescriptible attribute of the structure -nonlocality. Inapplicability of quasi-static approximation, that is effective medium model, wasn't demonstrated as well. We will start below from a demonstration of the fact that the effective medium model can't be applied for accurate description of the multilayered metal-dielectric metamaterials even in the deep subwavelength limit.

Dispersion Diagrams
First of all, we study electromagnetic properties of the lossless structure by means of dispersion diagrams and isofrequency contours. Losses will be introduced in Sec. 4.
For MMDN with two layers forming its period, using effective medium model and transfer matrix method we obtain two dispersion relations ω(k), where k = (k x , k y , 0) is the wave vector. TM polarization is assumed to take advantage of plasmonic behaviour. The first relation is given by effective medium model: and another one is given by transfer matrix method, assuming that k x is the wave number in the i-th layer: In this section graphic representation of the dispersion dependencies is chosen in the form of dispersion diagrams for waves travelling along y-direction, i.e. along the layers with k x = 0 and presented in Fig. 2. Isofrequency contours will be shown in Sec. 3.
General view of the structure is presented in Fig. 2.a. The structure consists of metal and dielectric layers and indefinite in all directions. Period of the structure is D = d 1 + d 2 = 62.5nm, where d 1 and d 2 are the layers thicknesses. Dielectric function for the dielectric layers is constant and equals to 4.6, while metal layers are described by dielectric function of the Drude form: with the plasma wavelength λ p = 2πc/ω p = 4D and the damping coefficient Γ. The multilayered metamaterial with 3 different ratios d 1 /d 2 of layers thicknesses (2/3, 1 and 3/2, respectively) and different values of D were chosen in order to demonstrate different behaviours of the structure. The case of the full period equal to 62.5 nm has been considered in Ref.
[29]. Here we also consider structures with periods equal to D/2, D/4, and D/8 and compare their behaviour with one of the full period. Astonishingly, the behaviour observed for the full period that is 10 times smaller than the wavelength remains qualitatively the same for structures having 2, 4 and 8 times smaller periods (i.e. 20, 40 and 80 times smaller than wavelength), as it is shown in Fig. 2. The convergence to an asymptote corresponding to surface plasmon-polariton (SPP) resonance becomes weaker as the branches reach asymptote at larger values of k y , but a topology of the branches remains the same. A pronounced tendency is such that the decrease of the period pushes one of dispersion curves closer to one predicted by the effective medium model for small k y . For large k y that convergence to effective medium model is weaker. Another dispersion branch tends to converge to a straight line corresponding to a volume plasmon with the decrease of the period. Volume plasmon is at the frequency where ε || = 0 and it is marked out with green thin line in Fig. 2. Positions of these two lines, effective medium model line and volume plasmon line, depend on the d 1 /d 2 ratio. That means that the convergence also depends also on the thicknesses ratio. In the case of equal thicknesses the convergence is faster than for different thicknesses. However, it is clear that even in the quasi-static limit correspondence between the effective medium model predictions and the actual dispersion properties of the metamaterial is not satisfactory, especially for high-k waves with relatively large k y which are mainly used in subwavelength imaging. We have also discovered from dispersion diagram analysis that the metamaterial in the case of thin metal (d 1 > d 2 ) possesses a fine structure shown in Fig. 2.d and supports four modes co-existing in the same frequency range.

Isofrequency Contours
Another powerful tool to analyze dispersion relation ω(k) is analysis of isofrequency contours, which are obtained if one fixes the frequency ω and plots k x (k y ) dependence. Isofrequency contours for three different thicknesses ratios of multilayered metamaterials at certain frequencies are presented in Figs. 3-5. You can find isofrequency contours for the whole optical frequency range in the video attachment of the present paper.
For the configuration corresponding to the case of thin metal (d 1 > d 2 ), which has the most unusual dispersion properties, isofrequency contours are presented in Fig. 3 for 8 different frequencies. Examination of the contours reveals that at certain frequencies they have very unusual topology. For D/λ = 0.066 (see Fig. 3.a), two hyperbolic isofrequency contours predicted by effective model are well reproduced by MMDN. Deviations from the effective medium model appear only at large k x close to the edges of the Brillouin zone, where the homogenization procedure is forbidden since the wavelength in the material becomes comparable to its period. For D/λ = 0.087 (see Fig. 3.b), where the frequency is closer to the SPP resonance, in addition to the hyperbolic contours MMDN demonstrates an elliptic-like contour right around the center of Brillouin zone. This region of the Brillouin zone is usually treated as one where the quasi-   static approximation is to be applicable since eigenwaves have small wave numbers, i.e. small spatial variation. However, inside layers fields vary dramatically and this leads to emergence of nonlocality. At slightly higher frequency D/λ = 0.089 (see Fig. 3.c), elliptical and hyperbolic contours observed at lower frequencies are transformed into a pair of crossing ellipses. In this case, the effective medium model predicts completely flat isofrequency contour which is perfect for realization of canalization regime for subwavelength imaging [5,8], but it is not in agreement with realistic dispersion properties of MMDN. At the frequency D/λ = 0.096 (see Fig. 3.d), an elliptical isofrequency contour predicted by effective medium model is accompanied by two hyperbolic contours, which disappear at higher frequencies, for example, at D/λ = 0.106 (see Fig. 3.e). However, they appear once again at slightly higher frequencies, e.g. D/λ = 0.106 (see Fig. 3.f). Presence of elliptic isofrequency contours is a clear sign of nonlocality in the MMDN. At D/λ = 0.108 (see Fig. 3.g) the topology of isofrequency contours is similar to ones at D/λ = 0.089 (the crossing ellipses), but at slightly higher frequency D/λ = 0.109 (see Fig. 3.h) the contours acquire even more unusual, pitcher-like form, which also has nothing to do with effective medium model predictions.
For the rest of thicknesses ratios when the thicknesses are equal or dielectric is thinner than metal isofrequency contours are presented in Figs. 4, 5. There are not new topological structures that are absent in the case of thin metal, however, hyperbolic-like and elliptical contours are still there. One should note flat isofrequency contours in Fig. 4c that appear at the frequency where ε ⊥ = ∞. It is known that flat isofrequency contours are in use for subwavelength image canalization [8,30]. Also, in the case of equal thicknessess plotted in Fig. 4 we see that isofrequency contours of effective medium model are in a quite good agreement with realistic ones.

Complex Band Structure
Dispersion relation given by Eq. 3 allows us to study complex modes of the structure immediately complex k y is considered. There are plenty of methods complex dispersion equation can be solved by, including a method based on the argument principle theorem [31], operator theory methods [32], and numerical FEM methods [33,34]. However, in order to find roots of Eq. 3 in the complex plane we employ a combination of change of sign method and Newton's one.
In Figs. 6, 7 dispersion diagrams of lossless (Γ → 0) and lossy structures showing complex modes are presented. We show positive real part of the wave vector with corresponding imaginary part. Due to the periodicity and unboundedness of the structure its complex eigenvalue spectrum contains an infinite number of modes. However, only modes having low imaginary parts of its wavenumber are of interest since they can travel significant distances. For example, modes having Im(k y ) < π/D are able to propagate to distances more than D/2π. We restrict imaginary part of k y by the value of 8π/D in our calculations.
In the absence of losses the structure supports the set of waves having the wavevectors ±k and ±k * . Thus, there are forward and backward waves along with their complex conjugations, and imaginary parts of the modes are symmetric in Figs. 6a, 7a, . It follows from the fact that if dielectric permittivities of layers are real then the right part of Eq.3 is a holomorphic function whose restriction to the real numbers is real-valued.
The most important three modes are marked in Figs. 6, 7 with different styles and colors while high-order complex modes are in gray color. Effective medium model results are presented in the figures as well (thin curves). Modes showed in Fig. 2 of Sec. 2 can be recognized by their zero imaginary part in regions of propagation. Where in Fig. 2 there were band gaps, here we have fully complex modes. At the SPP resonance frequency real part of the propagating mode goes to infinity, while imaginary part degenerates at this point what can be seen on the example of mode II in the case of d 1 < d 2 and on the mode I in the case of d 1 > d 2 . One can also observe an effect of decay switching in high-order modes: after a mode under- goes the SPP resonance its imaginary part jumps from one branch to another. In the case of d 1 < d 2 that is presented in Fig. 6a there is a very special mode III that is fully imaginary in the presented frequency range, and its real part is dispersionless and equals to zero. This feature is very unique since that mode is not sensitive to the resonance. Its presence is highly related to effective medium model as the latter has the resonance at another frequency and at higher frequencies converges to the entirely imaginary mode III while at low frequencies it converges to the mode II.
However, zero losses aren't realistic and unachievable in practice. Thus, we introduce attenuation in metal to study effect of losses in metal layers on electromagnetic properties of MMDN. Increasing losses smoothly we find that it changes step-like behaviour of the modes at the SPP resonance. That is, a specific form of these steps is determined by the value of losses in MMDN. With change of Γ modes go through a number of transformations [35]. Mechanism of modes transformation is exposed in Fig. 8. When losses tend to zero, at the resonance positive imaginary branches jump to their left neighbours with growth of the frequency, while negative branches jump to their right neighbours. Then, increase of the losses level forces branches to modify their topology and after a number of metamorphoses at the resonance both positive and negative imaginary branches jump to their right neighbours with growth of the frequency.
This way we come to realistic level of losses. Γ = 1.734 · 10 13 s −1 is supposed in our calcula- tions. Since now we analyze complex modes in the case of realistic losses comparing complex band structure with the lossless case to show how losses affect the behaviour of the structure.
In the presence of losses, as soon as the damping coefficient Γ in the dielectric function of metal ε 2 becomes positive non-zero and ε 2 becomes a complex number, the right part of Eq. 3 can be complex-valued with real k y . Consequently, the modes confluence is removed and the conjugated waves with ±k * are no longer supported by the structure. However, correspondence between positive and negative real and imaginary parts can be quite diverse because of the absence of symmetry in the imaginary parts of eigemodes (see Figs. 7b, 6b). In Figs. 7b, 6b we still have decay switching effect that is seen better in high degenerative modes. It's expected for modes we have observed in Figs. 7a, 6a: losses of a mode are grown up dramatically when it enters into a band gap, but for complex modes such behaviour is nonobvious. Also, as soon as one moves from the real axis to the complex plane every mode becomes complex one and it becomes hard to distinguish between propagating and evanescent modes. The case of d 1 > d 2 is shown in Fig. 6b. Just that very case corresponds to simultaneous presence of forward and backward modes in the same frequency range (D/λ : 0.105−0.108) for propagation along the layers. That is, the beam splitting phenomenon revealed in [20] remains even after introduction of realistic level of losses. Decay switching effects are observed there for high degenerative modes and the modes I, III while the mode II has narrow propagation domain with high losses in the rest of it.
Finally, it's noteworthy that effective medium model predicts resonance behaviour for imaginary part of k y . In reality, however, the structure is functioning in a regime with infinite number of modes which imaginary parts have step-like behaviour at the resonance.

Field Profiles
In order to obtain deeper understanding of the nature of the studied complex eigenmodes of the system, we have plotted profiles of the magnetic and electric fields of these modes. Fig. 9 shows the profiles of the 3 branches for the case of the thicker metal at the frequency D/λ = 0.0875. We first note that the magnetic field of the mode I averaged over each layer is equal to zero.   It means that the mode I is effectively longitudinal. Such modes are well known in plasma physics as Langmuir modes [36]. It has been also recently shown that such exotic modes may exist in the waveguides with metal cladding and core made of hyperbolic media [37]. We can see however, that the longitudinal nature of the mode holds only within the effective media approximation and the local magnetic field inside each layer does not vanish. Modes II and III are conventional transverse modes.
Another interesting feature that can be noted on the profiles is the complicated field distribution of the electric field inside the metal layer for mode I. Conventionally the modes existing in metal-dielectric structures are coupled surface plasmon resonances at the individual metaldielectric interfaces. Thus, there may exist no more than one extremum of the field distribution function in each layer. This condition however does not hold for the case of the complex modes. Indeed, if we consider the modes with real and imaginary part of the propagation constant, the expression for the transverse component of the wavevector in each layer reads: It is clear that when Im(k y ) is large enough, the k x gains large real part and the mode starts to propagate inside the structure. This leads to nonzero values of the transverse component of the Poynting vector inside the layers and to the complex shape of the field distribution inside the layers. To illustrate how the imaginary part of k y affects the shape of the field inside the layers we have plotted the profiles of the high order complex modes [see Fig. 10]. We can see that while we are switching to higher order modes, having the larger imaginary part of k y the field profile shape becomes more and more complicated. These modes can be regarded as a specific type of coupled waveguide modes. The waveguide mode condition can be roughly estimated as k i,x d i = πn, where i = 1, 2 and n -is an integer. Thus, higher switching between the complex modes can be regarder as switching between different waveguide modes in the structure.
We have also plotted the profiles of the three mode branches for the thinner metal case [see Fig. 11]. We can see that in this case we also have one longitudinal mode I and two transverse modes. In this case however the mode II has the complicated structure of the electric field as the one defined by the larger values of Im(k y ).

Conclusion
Using dispersion diagrams along with isofrequency contours we have analyzed different characteristic configuration of MMDN both in the absence and in the presence of losses. Study of real-valued eigenmodes has revealed the impossibility of its description in terms of the local effective medium model even in the deep subwavelength limit. With the introduction of realistic losses, complex band structure of MMDN was revealed. Generally, nonlocality in MMDN is slightly affected by losses. Thus, in real structures with losses variety of spatial dispersion effects is expected.