Three cases of discontinuous refractive index in metamaterial study

We investigate three cases of metamaterials presented in the literature displaying refractive index with one or more discontinuities along the frequency spectrum. We reproduce the numerical simulations of these metamaterials and compare our simulations to each reported case. For each case, we perform a geometrical investigation of each metamaterial’s refractive index by mean of a numerical simulation of a prism made of the reported metamaterials upon which is incident a plane electromagnetic wave. Such investigation allows us to infirm or confirm negative refraction at resonance frequency. Finally, we carry a numerical and theoretical investigation of this discontinuity and show that, as the refractive index crosses a discontinuity, while the topology of the effective wave has changed within the metamaterial, the dynamics of the phases remain unchanged at any time at the metamaterial's boundaries.


Metamaterials.
It was not until 2000 that Smith et al. 2 produced a composite medium displaying the double negative permittivity and permeability Veselago predicted. This material, by the fact that it was composite and displayed exotic characteristics, i.e. characteristics not found in nature, fitted the definition of metamaterials. This experimental realization caused an explosion of interest toward NIMs. A major figure in this initial impetus was John Pendry, who, a few months after Smith et al., suggested that a lens made of NIMs would beat the "natural" lenses resolution (i.e. no smaller than a wavelength) 3 , and in 2006 introduced, at the same time as, yet independently from, Ulf Leonhardt, the idea of a metamaterial-based invisibility cloak 4,5 . Such theoretical promises, quickly followed by experimental realizations [6][7][8] , reinforced the academic interest surrounding negative refraction. While the original design by Smith et al. was based on the so-called split ring resonator (SRR), a variety of designs appeared throughout the following years [9][10][11][12][13][14][15][16][17] . Most of the proposed NIMs, Zhang et al. 's fishnet included, are termed metasurfaces, presenting as films that can be laid on top of one another in order to create a layer of metamaterial of desired thickness 18 (see Fig. 1). Another often recurring specificity of NIMs is their anisotropy: most NIMs only display negative refraction when the plane wave incident upon them propagates www.nature.com/scientificreports/ normally to the surface and is polarized along one specific axis. These constrain the shapes and uses for which metamaterials can be used. 3D metamaterials do exist [19][20][21] , but have been much more challenging to develop. Another challenge is the limitations in transmission and frequency of operation of metamaterials. Most metamaterials developed so far are operable only in the infrared to microwave frequency range, as the development of negative refraction at optical frequencies meets serious technical challenges due to the dimensions constraints of the metamaterial's structure 22 . Metamaterials also display significant absorption around magnetic resonance; such absorption is often quantified using the Figure of Merit (FOM), calculated as the ratio of the real to imaginary part of the refractive index. Absorption is evidently an undesirable effect to any application of metamaterials: a metamaterial-based imagery device would reduce the quality of the resolved image, while an invisibility cloaking system, typically designed for stealth, could be detected by noticing the unusual absorption of the system's location.
Effective electromagnetic parameters retrieval. One way to retrieve the NIM's effective electromagnetic parameters, and in particular the refractive index, requires a knowledge of the complex transmission and reflection coefficients of the wave incident to the material. Regardless of whether transmission and reflection obtained by experiment or numerical simulation, one needs a way to convert these to the optical parameters. Such an algorithm has been developed. In 1970, Nicolson et al. 23 are the first to derive the set of equations relating transmission and reflection to refractive index and impedance, developing the well-known Nicolson-Ross-Weir (NRW) method. Chen et al. 24 in turn improved this algorithm, allowing the first boundary and thickness of the effective metamaterial (those being different from the physical boundaries and thickness of the metamaterial, as the electromagnetic wave is usually affected ahead and after the physical boundaries by the metamaterial) as well as the sign of the effective impedance and the branch of the real part of the refractive index to be determined. More recently, Hsieh and Wang 25 suggested another method applied to bianisotropic materials and based on material dispersion models. The most important assumption that the above algorithms rest upon is that of the material's homogeneity: the inhomogeneity within the material, mainly, its stratified composition and planar pattern, have to be designed within length scales that are very small compared to the effective wavelength of the wave propagating within the metamaterial (as an example, Ding et al. 's fishnet 14 is constituted of 30 μm-thick unit cells, displaying negative refractive index for electromagnetic waves of wavelength 0.1 mm, three times larger).
In their article Chen et al. point out how, around resonance frequency, the wavelength within the material approaches the material's length scale, thereby excluding the possibility of a homogeneity assumption and making the situation one that is outside of the range of application of their algorithm. This is true inasmuch as the wavelength within a material equals the product of the wavelength in free space and the refractive index of the material: the larger the refractive index, the smaller the wavelength within the material. Regardless, this homogeneity assumption has often been disregarded in other papers 11,[13][14][15][16][17] , who make the choice of presenting the results of the retrieval of effective electromagnetic parameters as accurate, even in the region of magnetic resonance, precisely where the wavelength of the radiation within the medium reduces to the same order of magnitude as the thickness of the material. The fact that the algorithm's hypothesis is not verified in the window of frequency where magnetic resonance occurs is tantalizing since this is precisely where the refractive index turns negative! Smith et al. 26 solved this problem by suggesting an algorithm to retrieve effective electromagnetic parameters in an inhomogeneous material. The 11,13-17 papers referenced above, however, do not refer to this method and insist on wrongly using Chen et al. 's algorithm. One more algorithmic challenge, to which Smith et al. has not brought an answer, is that of what we will hence-forth refer to as the m-branch index jump. In the above algorithms, the expression of the refractive index contains an inverse cosine, and is therefore dependent on an associated m-branch, as can be seen in Eq. (2): The determination of such m-branch index is usually done using physical requirements on the system or by a requirement that the real part of the refractive index be continuous. Some metamaterials 16,17,27 , however, have been reported to display a discontinuous refractive index, jumping from positive to negative, i.e. with the m-branch index jumping from one value to another. No consideration of this case seems to have been done when  www.nature.com/scientificreports/ developing the above algorithms, leaving the decision of the frequency at which the m-branch index jump occur up to the arbiter of the researchers, such as in 16 . Some papers 27 do acknowledge the limits of the algorithm and the impossibility to decide where the jump occurs.

Proposed solution.
In this article, we present three cases of metamaterials presented in the literature displaying refractive index with one or more discontinuities along the frequency spectrum. We reproduce the numerical simulations of these metamaterials and compare our simulations to each reported case. For each case, we perform a geometrical investigation of each metamaterial's refractive index by mean of a numerical simulation of a prism made of the reported metamaterials upon which is incident a plane electromagnetic wave. Such investigation allows us to infirm or confirm negative refraction at resonance frequency. Finally, we carry a numerical and theoretical investigation of this discontinuity and show that, as the refractive index crosses a discontinuity, while the topology of the effective wave has changed within the metamaterial, the dynamics of the phases remain unchanged at any time at the metamaterial's boundaries.

Three cases of discontinuous refractive index
In this section, we investigate what appears to be an anomaly of the NRW method: the discontinuity of the refractive index in some materials. We will first consider three designs presented in the literature, presenting their characteristic in details. We will then offer an alternative way to measure the refractive index of metamaterials, and suggest a possible explanation for the emergence of the discontinuity. We have identified three designs reported in the literature, which present discontinuous refractive index. Numerical simulation of these metamaterials is performed using the commercial software CST Microwave Studio 30 and results are compared to those of the respective reports, in order to confirm the faithfulness of our numerical reproduction of those metamaterials. In each case, a unit cell of metamaterial is built on the software, lying in the XY plane. Boundaries are set as magnetic (H t = 0) along the X-axis and electric (E t = 0) along the Y axis, with open boundaries along the Z axis, axis of propagation of the plane wave. A plane wave polarized along the Y axis is sent for normal incidence upon the unit cell, and the S-parameters are retrieved for analysis using the NRW method. Meshing is carefully monitored to always exceed 50,000 tetrahedrons.

Chan et al.'s arrays of upright split-ring pairs.
Chan et al. 17 proposed a metamaterial operating in near-infrared frequency, consisting in split-ring pairs arranged back to back. The split-ring pair is made of gold (Drude permittivity of ε = 9 − ω 2 p ω(ω+iŴ) ), with ω p = 1.37 × 10 16 Hz ) and Ŵ = 1.0027 × 10 14 Hz , while the encasing unit cell is in glass (n = 1.458; in our own simulations we will consider this to be fused quartz, as it fits the definition of glass and its refractive index matches that indicated by Chan et al. at near-infrared frequency 28 ).
A unit cell is represented in Fig  www.nature.com/scientificreports/ the effective parameter retrieval method, geometrical analysis of the incident rays and by analysis of the dispersion relation of their material. Nowhere in their report do they, however, mention the refractive index jump that occurs at 238 THz (from positive to negative value) and 297 THz (back from negative to positive value) (Fig. 3a).
It is also noteworthy that Chan et al. 's results are the only one among the three we will discuss here, where continuity is enforced-a change in the m-branch index of Eq. (2) at frequencies where the discontinuity can be observed need not take place and the retrieved refractive index can be made continuous across the spectrum of frequencies we're considering here. We chose to disregard this change in m-branch index in Fig. 3 to make the similarity of our numerical simulation to that of Chan et al. 's apparent. When the NRW algorithm is thoroughly applied and continuity of the refractive index is enforced, the refractive index becomes as represented in Fig. 4, and is positive for all frequency values of the considered spectrum.   Fig. 6, the application of the effective parameter retrieval procedure holds seemingly very different results: our application of the procedure leads to a negative refraction appearing only starting at 0.92 THz. A closer inspection reveals however that this is only due to the refractive index jump from positive to negative values, which occurs at an earlier frequency in Du et al. 's analysis. This discrepancy notwithstanding, the profile of the refractive index is very similar if one ignores the discontinuity and only considers the rate of change of the refractive index with respect to the frequency. This rate of change is positive from 0.6 to 0.76 THz, reaches a zero (refractive index local maximum) at 0.76 THz, then turns negative until reaching another zero (refractive index local minimum) at 0.81 THz, and then being positive again. This discrepancy in results highlights the importance of the phenomenon of refractive index discontinuity: Du et al. is able to report a frequency window of negative refraction about five times greater than that one would conclude from our own results. However, no acknowledgment was made of the refractive index discontinuity. Paul et al. 27 present a fishnet-like, cross-shaped metamaterial operating at low terahertz frequency. This design is similar to Du et al. 's excepted the presence of arms, which interestingly makes it the only "disconnected" fishnet considered in this article (see Fig. 7). The metallic layers are set to be copper, with Drude permittivity (the Drude characteristic of copper at terahertz frequencies are not specified in Paul et al. 's report; for our own simulation, we set those characteristics to ω p = 1.914 × 10 15 Hz for the plasma frequency and γ = 8.34 × 10 12 Hz for the collision frequency 29 ). Same as with Du et al. 's fishnet, the whole unit cell is permeated with the dielectric, chosen by Paul et al. to be BCB (permittivity of ε = 2.67 and loss parameter tanδ = 0.012 at frequency of interest).Paul et al. reports negative refraction in the frequency window from 0.96 to 1.17 THz. Our own simulation qualitatively agrees with those results, as can be seen in Fig. 8, albeit with a negative refraction starting at the refractive index discontinuity located at 0.82 THz. Paul et al. is the only one to actually acknowledge the discontinuity in their results, stating: "The dashed parts of the curves (Fig. 8a) mark the regions where the half-wavelength inside the medium is smaller than the lattice constant a z and the calculated effective parameters are not expected to be reliable" 27 .

Geometrical study of simulated electromagnetic dynamics as a test for negative refraction
A benchmark was used to investigate the negative refraction (or lack thereof) in metamaterials. This benchmark relies on the geometrical ray formalism of electromagnetic waves, using observations of the behavior of those waves as they are transmitted through objects of determined shapes (mainly, a prism), and comparing those observations to what would be expected should the material be homogeneous and have the claimed refractive index.
The numerical tools at our disposal for the simulation of plane waves as a way to observe the reaction of transmitted waves through the media. This brings a challenge to the analysis of the refractive index using this benchmark: diffraction patterns occur, that, while still allowing an interpretation of the sign of the refractive index, limit the precision of an exact computation of that index. Fortunately, in this analysis we are less concerned www.nature.com/scientificreports/ by the exact value of the index, than we are about its sign, and the frequency at which the jump from positive to negative refractive index occurs. The prism configuration is a benchmark that allows to geometrically determining the refractive index of the metamaterial (Fig. 9). A 2D prism is numerically designed using unit cells of the metamaterial one wishes to study, and a plane wave is sent onto a surface of the prism. The direction of the transmitted wave allows the determination of the sign of the refractive index of the metamaterial at that wavelength. Typically, the prism  www.nature.com/scientificreports/ is designed so that the interface of incidence of the plane wave is parallel to the wave front, and the prism is elongated enough to contain the principal front of the wave. A standard prism is depicted in Fig. 9. The unit cells (cyan) are stacked in a half-pyramid, creating a prism in the XZ plane. A plane wave (wave fronts depicted in black) polarized in the Y direction and propagating towards the Z direction is incident upon the lower interface of the prism, and emerges refracted from the upper interface  One readily notes a first challenge with this configuration: virtually all metamaterials are constrained by the dimensions of their unit cell, and the development of complex geometrical objects using those metamaterials must take these dimensions into account. For the examples considered in this paper, it is impossible to design a prism with an upper surface that can be considered nanoscopically or even microscopically planar, and that can therefore match the ideal upper interface; the best one can do is to design staircase-shaped prisms (grism) in dimensions such that a negligible amount of scattering occurring due to the local nonnuniformities of the upper surface. However, practically speaking this staircase shape does not lead to significant scattering if we constrain the slope by an increment of one unit cell's depth as we move laterally towards the positive X. An example demonstrating this can be found in Fig. 10, representing a prism with similar geometry as that of Fig. 9, made of 15(15+2) 2 homogeneous blocs of Teflon with dimensions 150 × 150 × 30 μm. An electromagnetic wave of frequency 1 THz, polarized along the Y-axis and propagating towards the positive Z, is sent upon the inferior interface of the prism. One clearly distinguishes the main refracted wave leaving the upper interface of the prism towards the right (as expected from a positive refraction), as well as secondary waves of lower intensities refracted right and left of the main wave, due to diffraction effects of the input port. Figure 11 shows a simulation of a plane wave incident upon a prism made of unit cells of Chan et al. 's metamaterial. The absolute value of the electric field is represented within and around the prism (the shape of which is explicit in black for the first frequency).

Analysis of Chan et al.'s metamaterial.
Comparing those results to those obtained by (mis-)use of the NRW algorithm, no negative refraction is observed at frequencies where expected from Fig. 3b. Instead, positive refraction is observed at frequencies 160-180 THz. At those frequencies, transmission drops to levels that prevent any analysis of the refracted plane wave. Between frequencies of 210-230 THz, strong reaction, confirmed to happen at those frequencies by numerical simulation of a single unit cell, causes the incident plane wave to "escape" laterally from the lower interface of the prism. These results confirm the necessity of enforcing continuity of the refractive index whenever possible. In the last two cases however, discontinuity has to occur, as we will show. Figure 12 shows a simulation of a plane wave incident upon a prism made of unit cells of Du et al. 's metamaterial. The absolute value of the electric field is represented within and around the prism.

Analysis of Du et al.'s fishnet three-dimensional metamaterial.
Positive refraction is observed from 0.6 to 0.72 THz, after which a drop in transmission is observed (and confirmed in the numerical results presented in Fig. 6, with a transmission's amplitude under 0.60 in the span 0.73-0.82 THz and a minimum of 0.10 reached at frequency 0.78 THz), rendering any observation impossible. Negative refraction is observed at 0.86 THz transmission rebounds. Enhancement of the electric field is done by rescaling of the color gradient. However, because of the high loss of the material only the part of the plane wave that crosses a small thickness of prism emerges from the upper interface, skewing the results (see Fig. 13 for an illustration of the effect; the white star is the desired center of the wave emerging from the upper interface, and the black star the observed center). Thus, this analysis does not allow for an exact determination of the frequency at which the refractive index passes from positive to negative, but allows to determine a region of uncertainty. In particular, it points out at an error in our application of the NRW algorithm, as our computation of the refractive index predicted a negative refraction only starting frequency 0.92 THz. Details about the origin of this error will be covered later in this article. Positive refraction is observed from 0.6 to 0.74 THz. Starting 0.86 THz, a clear negative refraction is observed. Between 0.74 and 0.86 THz, an interesting phenomenon occurs: the emerging positively refracted wave front gradually diminishes in intensity while a negatively refracted wave front gradually increases in intensity, with a seemingly equal wave intensity observed at frequencies 0.80-0.82 THz, where the plane wave seems to split at the upper interface.
Regardless of the ambiguity observed between frequencies 0.74-0.86 THz, our observation allows us to extend the reported frequency spectrum over which negative refraction occurs (the grey zone in Fig. 8a) from starting at 0.96 THz to starting at 0.86 THz.  Figure 15 shows, side by side, the numerically retrieved transmission phase of the metamaterials, as well as the numerically retrieved refractive index with, in the case of Chan et al., misuse of the NRW algorithm (i.e. non enforcement of the continuity requirement). One readily notices that refractive index discontinuities occur at the exact same frequencies as a "modulo jump" of  In other words, the real part of the refractive index possesses a term that is directly proportional to the transmission phase. One notices a correspondence between this result and Fig. 15: the jumps in transmission phase and refractive index are opposite in direction, one going up when the other goes down, justified by the negative sign of the arg[S21] term in Eq. (5). Furthermore, these phase jumps correspond exactly to a m-branch jump, since incrementing the m-branch index by one is equivalent to dropping the transmission phase by 360 degrees (2π). This is therefore proof that close monitoring of the m-branch index is able to solve mistaken discontinuities in the application of the NRW algorithm, such as with Chan et al. 's.
One would be tempted to formulate a rule to always accompany a transmission phase jump by a corresponding change in m-branch index, thereby always enforcing continuity. This, however, works only in the case where enforced continuity is validated by numerical simulations such as tests carried by prism refraction. In the cases  One notices the term associated with spatial propagation has become negative (bearing in mind that d ≪ λ meta ), as expected from the ip in wave vector k direction that occurs when passing from positive refraction to negative refraction (i.e. the wave vector k goes in the same direction as the Poynting vector in regular materials, but in the opposite direction in metamaterials). Equation (7) reveals that the phases at any boundary of a stack of metamaterial unit cells remain unchanged after a refractive index jump, specifically being incremented by a multiple of 2π: considering a stack of thickness md with m a signed integer (including 0) and d the thickness of a single unit cell, the phase at z = md passes, as the refractive index jumps to a lower m-branch index, from mdk − ωt to mdk − 2 mπ − ωt. However, the sine wave within the metamaterial has been flipped," as shown in Fig. 16. The angular frequency ω changes neither in value nor in sign, so that the wave internal to the metamaterial has now a negative phase velocity.
These considerations show that, while the topology of the effective wave has changed within the metamaterial, the dynamics of the phases remain unchanged at any time at the metamaterial's boundaries. Bearing in mind that the effective wave described by the effective parameters resolved by application of the NRW algorithm is just a convenient model of the metamaterial's properties, and not and actual description of the electromagnetic dynamic within, so that no actual "wave flip" occurs, the refractive index discontinuity is therefore entirely compatible with upon which a plane wave is incident.

Conclusion
In this article, we have analyzed three instances of metamaterials reported in the literature who display a discontinuity on the frequency spectrum. We performed numerical simulations of metamaterial unit cells with similar parameters as those described in those reports to ascertain a faithful reproduction of the report's research, and match the results of our numerical simulations with those of those reports. We then performed a numerical simulation of a prism constituted of the aforementioned metamaterials, upon which was incident a plane wave of frequency of resonance of the metamaterial. We were thus able to identify one metamaterial that did not display negative refraction, and two were. In particular, we were able to identify that Paul et al. 's metamaterial displays a beam-splitting phenomenon that warrants further research.
We identified further the cause leading to the emergence of the refractive index discontinuity, and developed a mathematical formalism to justify the connection between this discontinuity and the 2π modulo correction of the transmission phase. We finally demonstrated that while this discontinuity leads to a change of the behavior of the effective electromagnetic wave within the metamaterial, the behavior of the electromagnetic wave external to the metamaterial remains unchanged; as this discontinuity does not affect the phase of the wave at the metamaterial's boundaries. www.nature.com/scientificreports/

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.