Klein tunneling near the Dirac points in metal-dielectric multilayer metamaterials

The Klein tunneling of optical waves near the Dirac points in the metal-dielectric multilayer metamaterials is theoretically investigated and demonstrated through the coupled-mode theory under the tight-binding approximation and the rigorous band structure analysis based on the transfer-matrix method. The optical analogue of Klein tunneling for the relativistic electrons passing across a potential barrier is revealed by the iso-frequency contour analysis and numerical simulation to describe the optical beam propagation and refraction across the interface of two metal-dielectric multilayer metamaterial stacks. The transmission and reflection spectra of the Klein tunneling of optical waves are also explained by the coupled mode theory.

Klein tunneling represents the phenomenon that the relativistic fermions can pass across a high potential barrier without the exponential damping, which is a particular property that arising from the existence of negative-energy solutions of the Dirac equation [1][2][3][4][5] . Although, the Klein tunneling for charge carriers acting as massless Dirac fermions in graphene heterojunctions is demonstrated [6][7][8] , the Klein tunneling for relativistic electrons is very difficult to be directly observed in experiments due to the required experimental condition of extremely high energy. Since Klein tunneling can be described as a wave phenomenon, optical analogues of Klein tunneling for light beam propagation in waveguide lattices and photonic crystals with conical dispersion are studied in numerical simulations [9][10][11][12] . Recently, metal-dielectric multilayer metamaterials with unique dispersion relation are utilized for realizing many intriguing optical phenomena such as negative refraction 13,14 , sub-wavelength focusing 15,16 , diffraction-free propagation [17][18][19] , and anomalous indefinite cavities 20 . The strong spatial dispersion in metal-dielectric multilayer metamaterials will induce the optical nonlocality, which is connected to the coherent coupling of surface plasmon polariton (SPP) eigenmodes propagating along the interfaces of metal and dielectric layers 21,22 . It is shown that the degeneracy of the symmetric and anti-symmetric SPP modes in metal-dielectric multilayers forms the Dirac point at the centre of the Brillouin zone, where an effectively zero "optical mass" is realized 23,24 . The existence of the Dirac points in metal-dielectric multilayer metamaterials not only enable important research topics about the epsilon-near-zero metamaterials 23,24 and the nonlocal effective medium theory [25][26][27][28][29] , but also open new opportunities in connecting the classical electrodynamics system and the relativistic quantum mechanics system 30,31 for exploring relativistic phenomena for massless Dirac fermions in condensed matter physics such as Zitterbewegung effect 32,33 , etc. Furthermore, it is worth mentioning that except for the metal-dielectric multilayer metamaterials, the Dirac point can also be found in a simple lattice of graphene sheets embedded in a dielectric host medium, which can raise the topological phase transition in the system 34 .
In this work, the Klein tunneling of optical waves near the Dirac points in the lossless metal-dielectric multilayer metamaterials is investigated and demonstrated through theoretical analysis and numerical simulation. Compared with the binary dielectric superlattices that previously applied to study the Klein tunneling for optical waves 9, 10 , the lossless metal-dielectric multilayer metamaterials not only have a simpler structure, which can be tuned to achieve the Klein tunneling easily, but also possess a special band structure associated with the two SPP eigenmodes, which can be degenerated to form the Dirac point, resulting in a better analogue to the Klein tunneling in graphene. According to the specified coupled-mode theory under the tight-binding approximation, the connection between the metal-dielectric multilayer metamaterials and the relativistic fermions system described by the Dirac equation is revealed. The band structures of metal-dielectric multilayer metamaterials are analyzed based on the transfer-matrix method to reveal both the frequency regions of Klein tunneling around the Dirac points and the Klein tunneling zone in the wave vector space. The optical analogue of Klein tunneling for the relativistic electrons passing across a potential barrier is further illustrated by the iso-frequency contour analysis to describe the optical wave propagation and refraction across the interface between two metal-dielectric multilayer metamaterial stacks, which is coincident with the optical beam tunneling visualized from numerical simulation. Furthermore, the transmission and reflection spectra of the Klein tunneling of optical waves across the interface of two multilayer metamaterial stacks are also explained based on the coupled-mode theory.

Results
Coupled-mode theory and Dirac equation. As displayed in Fig. 1(a), two metal-dielectric multilayer metamaterial stacks with an interface located at x = 0 are considered. The multilayer metamaterial stack along the negative x-direction in Region 1 is composited of infinite alternating layers of gold (Au) and fused silica (SiO 2 ), while the multilayer metamaterial stack along the positive x-direction in Region 2 is made of infinite alternating layers of gold (Au) and alumina (Al 2 O 3 ). The light propagating inside the metamaterials in the x-z plane is TM (transverse magnetic) polarized. The thickness of the Au layer is a m = 10 nm, while the thicknesses of the SiO 2 layer and the Al 2 O 3 layer are the same as a d = 400 nm. The permittivity of Au is describe by the Drude model as with the offset constant ε ∞ = 9 and the plasma frequency ω p = 13.8 × 10 12 rad/s 35 . The damping factor γ in the Drude model is set to be zero in the current analysis in order to emphasize the SPP eigenmodes that associated with the Klein tunneling. In addition, the permittivities of SiO 2 and Al 2 O 3 are simplify set as ε = . × . , respectively. Note that in order to distinguish the same physical quantities in either the Au-SiO 2 multilayer or the Au-Al 2 O 3 multilayer, the superscripts of (1) and (2) are applied.
As depicted in Fig. 1(b), each metal-dielectric multilayer stack can be regarded as the combination of infinite coupled dielectric-metal-dielectric waveguide array, where a single waveguide with two eigenmodes can be numbered by the index of n, while the two eigenmodes supported in the nth waveguide are numbered as the 2n mode and the 2n + 1 mode, respectively 38 . Therefore, the band structure of each metal-dielectric multilayer stack can be described by the coupled-mode theory 39 . Under the tight-binding approximation, the coupled-mode theory for each of the two metal-dielectric multilayer stacks can be expressed as where β 0 presents the propagation constant of optical wave, δ presents the propagation constant difference of mode, and κ presents the mode coupling strength [Appendix 1]. Moreover, the wave functions A 2n and A 2n+1 individually denote the 2n mode and the 2n + 1 mode in the n th waveguide, while the wave functions A 2n−1 and A 2n+2 individually denote the 2(n − 1) + 1 mode in the (n − 1)th waveguide and the 2(n + 1) mode in the (n + 1) th waveguide. According to the Bloch's theorem, Eqs (1) and (2) can be solved by setting the wave functions A 2n−1 , A 2n , A 2n+1 , and A 2n+2 as with the periodic phase shift φ = k x (a m + a d ) and the amplitudes A and B respectively for the two modes in a single waveguide. With respect to the wave functions in Eq. (3), the band structure of the metal-dielectric multilayer stacks can be obtained from Eqs (1) and (2) [Appendix 1] as In general, Eq. (4) indicates that the band structure possesses both the β + -branch and the β − -branch with respect to the frequency ω. The relation between the propagation constants β ± and k x at a certain frequency ω based on Eq. (4) represents the iso-frequency contour (IFC). Since the tight-binding approximation only considers the coupling among the adjacent waveguides, the band structure predicted in Eq. (4) is only accurate around the Brillouin zone centre compared with the rigorous results based on the transfer-matrix method (see the following analysis). Nevertheless, Eqs (1) and (2) can still reveal the in-depth mechanism about the propagation of the TM polarized light in the metal-dielectric multilayer stacks. By introducing the relations of / under the assumption of the slow mode amplitude variation between the adjacent waveguides, where a = a m + a d is the thickness of multilayer unit cell, Eqs (1) and (2) can be rewritten as are the Pauli matrices. By comparing with the one-dimensional Dirac equa- it is shown that the coupled-mode theory in Eq. (7) has the same format as the Dirac equation with the coefficient mapping relations of κ ↔ c, δ ↔ mc 2 /ħ, and −β 0 ↔ V(x)/ħ, and the variable mapping relations of z ↔ t and ξ ↔ x. Furthermore, the positive-and negative-energy solutions of the Dirac equation can be mapped into the β + -and β − -branches of the band structures in the metal-dielectric multilayer stack, respectively. Therefore, the Klein tunneling for the relativistic electrons passing across a potential barrier in the relativistic quantum mechanics can be mimicked by its optical analogue where the TM polarized light waves propagate and refract across the interface of two metal-dielectric multilayer metamaterial stacks.
Band structure analysis for Klein tunneling. In order to precisely predict the Klein tunneling of optical waves across the interface of two metal-dielectric multilayer metamaterial stacks, the rigorous band structure based on the transfer-matrix method is considered to describe the propagation of the TM polarized light in each metal-dielectric multilayer stack as Generally, the rigorous band structure in Eq. (9) implies the relations among the wave vectors k x , k z , and the frequency ω, containing the information of both the IFC at a certain frequency ω and the dispersion relation at a certain wave vector k x (especially at the Brillouin zone centre k x = 0). Similar to the coupled-mode theory, the IFCs obtained from Eq. (9) for each metal-dielectric multilayer stack also possess two branches, in which the symmetric mode with wave vector β s and the anti-symmetric mode with wave vector β a can be defined at the Brillouin zone centre. Furthermore, the degeneracy of the symmetric mode and the anti-symmetric mode forms the Dirac points, leading to giant optical nonlocality in the metal-dielectric multilayer stacks. Figure 2(a) displays the dispersion relation between the normalized wave vector k z /k p and the normalized frequency ω/ω p at k x = 0 for two metal-dielectric multilayer metamaterial stacks, where the plasma wave vector k p = ω p /c is associated with the plasma frequency ω p and the speed of light c in free space. The dispersion relations for the Au-SiO 2 and Au-Al 2 O 3 multilayer stacks are denoted by the blue curves and the red curves, respectively, including both the symmetric modes β s (1) and β s (2) and the anti-symmetric modes β a (1) and β a (2) . Clearly, the degeneracy of the symmetric and anti-symmetric modes leads to the emergence of the Dirac points at the frequencies ω D (1) and ω D (2) for the Au-SiO 2 and Au-Al 2 O 3 multilayer stacks, respectively. According to the theory describing the Klein tunneling in the relativistic quantum mechanics, the frequency band for the Klein tunneling of optical wave to occur is from ω L = 0.0788ω p (173.098 THz) to ω H = 0.120ω p (262.559 THz) for the current two multilayer stacks, where ω L is determined by the intersection of the dispersion curves for the symmetric mode β s (1) and the anti-symmetric mode β a (2) , and ω H is located at the intersection of the dispersion curves for the anti-symmetric mode β a (1) and the symmetric mode β s (2) . Furthermore, as shown in Fig. 2(a) the frequency band of Klein tunneling is divided into three frequency ranges as I, II, and III, with respect to the locations of two Dirac points. Figure 2(b) further shows the typical IFCs at the frequency of ω = 0.095ω p (208.652 THz) in the frequency range II for the Au-SiO 2 (solid blue curves) and Au-Al 2 O 3 (solid red curves) multilayer stacks. For the Au-SiO 2 multilayer stack, it is clear that the IFCs possess two different branches with the anti-symmetric mode β a (1) in the lower branch and the symmetric mode β s (1) in the upper branch at the Brillouin zone centre. The similar IFC profiles are also observed for the Au-Al 2 O 3 multilayer stack, but with the β a (2) mode in the upper branch and the β s (2) mode in the lower branch. This is because the selected frequency ω = 0.095ω p is located between the two Dirac point frequencies ω D (1) and ω D (2) for the two multilayer stacks and there is a band inversion across the Dirac points, which is coincident with the dispersion relations shown in Fig. 2(a). Besides, the IFCs obtained from the coupled-mode theory (dashed curves) are also plotted for comparison, indicating that the coupled-mode theory only works well around the Brillouin zone centre. Note that both the energy and momentum are conserved in the Klein tunneling for the relativistic electrons passing across a potential barrier. For the optical analogue of Klein tunneling working at a specific frequency ω, the wave vector k z is conserved as the light waves refracted at the interface of two multilayer stacks. Therefore, as indicated in Fig. 2(b), the Klein tunneling zone in wave vector space can be directly obtained from the overlapped IFCs of two multilayer stacks sharing the same wave vector k z .  -branch, as well by the coupled-mode theory with the overlap of the β + (1) -branch and the β − (2) -branch. It is worth emphasizing the Klein tunneling zones predicted from two approaches are different from each other, but they are coincident with each other at the Brillouin zone centre.
The IFCs obtained from the transfer-matrix method are used next to analyze the refraction processes at the interface of two multilayer stacks in the Klein tunneling of optical waves. For instance, Fig. 3 shows the IFCs at three different frequencies in the frequency range I (ω = 0.083ω p ), II (ω = 0.095ω p ), and III (ω = 0.112ω p ) based on the transfer-matrix method. As depicted in Fig. 3(a), The IFCs at the frequency of ω = 0.083ω p (182.296 THz) in the frequency range I indicate that an incident TM polarized light beam from air with 17.5° angle of incidence will excite two propagation modes within the Au-SiO 2 multilayer stack, where the first mode belongs to the β s   (1) and (2). The frequency band ω ∈ [ω L , ω H ] for the Klein tunneling of optical waves is depicted as three ranges as I, II, and III. (b) The IFCs at the frequency ω = 0.095ω p (208.652 THz) in the frequency range II for the Au-SiO 2 multilayer stack (blue curves) and the Au-Al 2 O 3 multilayer stack (red curves) according to the transfer-matrix method (solid curves) and the coupledmode theory (dashed curves). The Klein tunneling zones in wave vector space based on the transfer-matrix method (TMM tunneling zone) and the coupled-mode theory (CMT tunneling zone) are also indicated.
the Au-Al 2 O 3 multilayer stack and excite the mode with the beam propagation direction of 36.5°. In order to further visualize the Klein tunneling of optical waves through two multilayer stacks, Fig. 3(d) displays the numerically simulated optical beam propagation and refraction processes inside the multilayer stacks. The simulation shows that an incident TM polarized Gaussian beam from air with 17.5° angle of incidence (k x = 0.025 k p ) is split into two propagating beams in the Au-SiO 2 multilayer stack with one beam in the beam propagation direction of 14.8° tunneling into the Au-Al 2 O 3 multilayer stack and the refraction angle is 36.5°. The similar optical mode propagation and tunneling processes for the Klein tunneling at the other two frequencies of ω = 0.095ω p (208.652 THz) and ω = 0.112ω p (245.99 THz) can also be found in Fig. 3(b,c,e,f), respectively. In addition, it is shown that the Klein tunneling of optical waves can take place between optical modes with different symmetries, depending on the operation frequency. Generally, in the frequency range I, the Klein tunneling arises from the mode tunneling between the symmetric β s -branch. For comparison, the propagation of the TM polarized light outside the Klein tunneling ranges in either the wave vector space or the frequency domain is also studied in Fig. 4. The first case displayed in Fig. 4(a,c) demonstrates the beam propagation at the frequency of ω = 0.095ω p (208.652 THz) which is located in the frequency range II for the Klein tunneling, but with the wave vector k x = 0.1 k p out of the tunneling zone. The IFCs in Fig. 4(a) indicate that an incident TM polarized light from SiO 2 with 46.5° angle of incidence (wave vector k x = 0.1 k p ) excites the mode in the beam propagation direction of 12.5° in the β s (1) -branch of the Au-SiO 2 multilayer stack. However, since no mode exists in the β s (2) -branch of the Au-Al 2 O 3 multilayer stack to match the wave vector k z , the Klein tunneling cannot occur. Correspondingly, Fig. 4(c) shows that the optical beam will not refract across the interface of two multilayer stacks based on the simulation, which is coincident with the IFC analysis. However, it is worth mentioning that there is still some optical energy evanescently couples cross the interface but it gets decayed quickly. The second case shown in Fig. 4(b,d) studies the beam propagation at the frequency of ω = 0.125ω p (274.542 THz) which is above the frequency range III for the Klein tunneling. From the IFCs in Fig. 4(b), it is shown that an incident TM polarized light from air with 16.3° angle of incidence (wave vector k x = 0.035 k p ) excites two modes in the Au-SiO 2 multilayer stack, with the beam propagation direction of 25.8° in the β s  -branch of the Au-Al 2 O 3 multilayer stack with the equal wave vector k z , the Klein tunneling cannot occur. The simulation in Fig. 4(d) also clearly shows that the propagation mode with the direction of 11.2° in the Au-SiO 2 multilayer stack is completely reflected back at the interface of two multilayer stacks.

Transmission and reflection of Klein tunneling. Analogous to the Dirac equation describing the Klein
tunneling for the relativistic electrons, the coupled-mode theory reveals the propagation and refraction processes of the TM polarized light waves tunneling across the interface of two multilayer stacks. Moreover, the coupled-mode theory can also be used to determine the transmission and reflection properties of the Klein tunneling of optical waves. With respect to the two multilayer stacks, the transmission coefficient t and the reflection coefficient r can be derived based on the coupled mode-theory as x (1) with the parameters (1) and , the mode propagation constants β 0 (1) and β 0 (2) , the propagation constant differences δ (1) and δ (2) , and the coupling strengths κ (1) and κ (2) , for the Au-SiO 2 and Au-Al 2 O 3 multilayer stacks, respectively [Appendix 2]. Furthermore, compared to the relativistic quantum mechanics theory for Klein tunneling, the intensity current for the TM polarized light in either of the two multilayer stacks can be defined as Therefore, the transmission and reflection of the Klein tunneling read as According to Eqs (13) and (14), Fig. 5 plots the transmission and reflection spectra as a function of wave vector k x , together with the corresponding IFCs based on the couple-mode theory at three different frequencies considered in Fig. 3. The IFCs in Fig. 5(a-c) are obtained from the coupled-mode theory with both the β + and β − branches denoted for each multilayer stack, in which the two branches involved in the Klein tunneling are denoted by solid curves, while other branches are denoted by dashed curves. The boundaries of wave vector k x for the Klein tunneling zone are also marked in the IFC figures. As shown in Fig. 5(d-f), the calculated transmission increases from zero to the maximum and then goes back to zero once the wave vector k x is out of the Klein tunneling zone, where the incident light is totally reflected by the interface of two multilayer stacks. Correspondingly, the reflection just behaves in an opposite way. The transmission and reflection spectra match with the Klein tunneling zone predicted from the IFC figures. In addition, the results of transmission and reflection spectra also imply the relation of T + R = 1 under the lossless condition due to the energy conservation.

Discussion
The Klein tunneling of optical waves across the interface of two lossless metal-dielectric multilayer metamaterial stacks is demonstrated through the theoretical analysis based on the coupled-mode theory and the transfer-matrix method, together with numerical simulation. The connection between the metal-dielectric multilayer metamaterial stack modeled with the coupled-mode theory and the relativistic fermions system described by the Dirac equation is revealed. The Klein tunneling of optical waves is clearly illustrated by the iso-frequency contour analysis and the numerical simulation to visualize the beam propagation and refraction processes across the interface of multilayer stacks, together with the calculated transmission and reflection spectra. The demonstration of the Klein tunneling of optical waves in metal-dielectric multilayer metamaterials near the Dirac points due to the coupling between surface plasmon polariton eigenmodes will create new frontiers in exploring the Dirac-cone physics associated with relativistic particle behaviors in condensed matter systems.

Method
Coupled-mode theory. In general, the coupled-mode theory can be expressed as (15) k k z i ( ) which describes the relation between the mode ∼ ν A with wave vector k ν along the z-direction and all other modes ∼ µ A with wave vector k μ along the same direction through the coupling strength κ ν µ . Mathematically, Eq. (15) can be rewritten as (16) by introducing the inverse Fourier transform as ). Regarding the coupled dielectric-metal-dielectric waveguide array depicted in Fig. 1(b), the coupled-mode theory in Eq. (16) reads standing for the mode coupling within the same waveguides, Eqs (17) and (18) can be written as with respect to the propagation constant β 0 = (β 2 + β 1 )/2 and the propagation constant difference δ = (β 2 − β 1 )/2. By substituting the symbolic solutions in Eq. (3), the band structure (i.e., the eigenvalues) can be obtained as x 0 0 in in 2 2 with σ = κ out /κ in based on Eqs (19) and (20). Note that the band structure in Eq. (21) must be coincident with the results of the results of the rigorous solution in Eq. (9) at the Brillouin zone centre, thus following the simple definitions that β − (0, β 0 ) = β 1 = min(β s , β a ) and β + (0, β 0 ) = β 2 = max(β s , β a ), there will be σ = κ out /κ in = −1 and κ in > 0 for Eq. (21). Simply speaking, since the coupled dielectric-metal-dielectric waveguide array is constructed by infinite identical single waveguide, the coupling strength κ in and the coupling strength κ out can be treated to have the same in magnitude but with a fixed phase shift. Therefore, by defining κ in = −κ out = κ > 0, the coupled-mode theory of Eqs (19) and (20) will be deduced into the simple format of Eqs (1) and (2) (2) ( 2) (2) (2) 2 ( 2) for the TM polarized light. Meanwhile, at the interface of the two different multilayer stacks, i.e., the interface between the −1th waveguide and the 0th waveguide (note that = ⌊ ⌋ n x a / ), the amplitudes of the TM polarized light should also be conserved, which means (1) (1) (2) ( 2) 0 (2) Based on Eqs (24) and (25), the transmission coefficient and the reflection coefficient can be obtained as in Eqs (10) and (11), respectively.