Giant optical nonlocality near the Dirac point in metal-dielectric multilayer metamaterials

The giant optical nonlocality near the Dirac point in lossless metal-dielectric multilayer metamaterials is revealed and investigated through the analysis of the band structure of the multilayer stack in the three-dimensional omega-k space, according to the transfer-matrix method with the optical nonlocal effect. The position of the Dirac point is analytically located in the omega-k space. It is revealed that the emergence of the Dirac point is due to the degeneracy of the symmetric and the asymmetric eigenmodes of the coupled surface plasmon polaritons. The optical nonlocality induced epsilon-near-zero frequency shift for the multilayer stack compared to the effective medium is studied. Furthermore, the giant optical nonlocality around the Dirac point is explored with the iso-frequency contour analysis, while the beam splitting phenomenon at the Dirac point due to the optical nonlocal effect is also demonstrated.


I. INTRODUCTION
Dirac point, conical singularity first discovered in the band structure of graphene [1], is the source of many remarkable wave transport properties in both electronic [2][3][4] and photonic system [5][6][7][8][9], such as the long-range flavor current splitting induced by the giant nonlocality in graphene [10], conical diffraction [11], Bragg-like reflection [12], quantum-like Goos-Hänchen shifts [12], and photon Zitterbewegung [13,14] in photonic crystals, because of the linear dispersion relation in the neighborhood of the Dirac point.In addition, rigorous analysis based on the tight-binding approximation and group theory clearly reveals the necessary conditions [15,16] for the creation of the Dirac point in the photonic system, including the optical metamaterials, a kind of artificial composites with periodic subwavelength metaatoms and tunable electric permittivity and magnetic permeability.Particularly, analogous to the effectively massless electron induced by the linear dispersion near the Dirac point of the electronic system, a Dirac point located at the center of the Brillouin zone (BZ) in the photonic system and optical metamaterials implies an effectively zero "optical mass", i.e., zero electric permittivity [17,18], leading significant analogies between the propagation of light in such media and charge transport in graphene [19].
On the other hand, epsilon-near-zero (ENZ) metamaterials with exotic electromagnetic properties have been studied in both theory [20][21][22] and experiment [23][24][25][26] recently.It is generally thought that the ENZ regime in the domain of the optical frequency range guarantees the local response described by the effective medium theory (EMT), since the vanishing effective permittivity implies a long effective wavelength resulting in a negligible optical size of the period.However, associated with periodic meta-atom structures in metamaterials, the variation of the electromagnetic field at the scale of a single meta-atom will result in the optical nonlocality, which is studied in both nanolayered [27][28][29][30] and nanorod [31,32] metamaterials, and extended to transformation optics [33].Especially, the optical nonlocality is theoretically analyzed in the metal-dielectric multilayer metamaterials since the geometrical simplicity allows the exact analytical calculation [30].It is worth revealing the mechanism of the optical nonlocality near the Dirac point in the ENZ regime.
In this work, we demonstrate the giant optical nonlocality near the Dirac point in lossless metal-dielectric multilayer metamaterials.Based on the transfer-matrix method for onedimensional layered photonic crystals, the band structure of the multilayer stack is fully explored and illustrated in the three-dimensional ω-k space.The exact location of the Dirac point in the band structure is located by rigorous algebraic analysis.It is revealed that the degeneracy of the symmetric and the asymmetric eigenmodes of the coupled surface plasmon polaritons (SPPs) induces the emergence of the Dirac point in the band structure of the multilayer stack, based on the dispersion relation analysis and the mode analysis around the Dirac point.Moreover, the locations of the Dirac point and the ENZ response in the ω-k space for the multilayer stack with the optical nonlocality and the effective medium are also studied.The giant optical nonlocality near the Dirac point is revealed by analyzing the isofrequency contours (IFCs) in the k-space with respect to the EMT and the transfer-matrix method.Furthermore, the beam splitting phenomenon induced by the giant nonlocality at the Dirac point, which is the optical analogue to the long-range flavor current splitting at the Dirac point in graphene [10], is also demonstrated by numerical simulations.

II. BAND STRUCTURES AND DIRAC POINTS OF MULTILAYER STACK
Consider a multilayer stack composite of alternating layers of gold (Au) and alumina (Al 2 O 3 ), with the thickness of each layer as d 1 = 15 nm and d 2 = 35 nm, respectively, as depicted in Fig. 1(a).The permittivity of Au is described by the lossless Drude model p /ω 2 , with the permittivity constant ε ∞ = 5.7, and the plasma frequency ω p = 1.3666 × 10 16 rad/s.Here the material loss of Au is neglected, since it will suppress the optical nonlocality of the multilayer stack.The permittivity of Al 2 O 3 is ε 2 = 2.4.
The multilayer stack can be regarded as a homogenous effective medium with the effective permittivity components described by the EMT as It is worth mentioning that the EMT, Eqs. ( 1) and (2), only take into account the filling ratios of different materials in the multilayer stack, i.e., the filling ratio of Au as and the filling ratio of Al 2 O 3 as rather than the actual thickness of each layer.Hence, EMT is an approximated theory that works well when the period the multilayer stack is much smaller than the wavelength of the propagated electromagnetic wave within the metamaterial, in the condition of the long wavelength limit.
More precisely, the stack can be regarded as a layered photonic crystal with the period of d along the x-direction, and the period of infinity along the y-and z-direction.Regarding a TM-polarized electromagnetic wave (with non-vanishing components of E x , E y , and H z ) propagating in the x-y plane, the band structure of the multilayer stack reads according to the transfer-matrix method [27].Here k where k = ω/c is the wave vector in free space.The band structure of the multilayer stack calculated from Eq. ( 3) is displayed in Figs.1(b) and 1(c) in the first BZ with respect to the period d along the x-direction.It is clear that the band structure contains three separated bands [Fig.1(b)], while the band 1 and the band 2 are connected with each other only at two points [Fig.1(c)], i.e., the Dirac points that will be focused on here.
Regarding a simple dielectric-metal-dielectric, or metal-dielectric-metal multilayer structure, it is known that the symmetric and the asymmetric SPP modes are the two fundamental propagating modes in the structure.Meanwhile, in the band structure of the multilayer stack, the emergence of the two connecting points, i.e., the Dirac points, implies the mode degeneracy of the symmetric and the asymmetric eigenmodes.Therefore, in order to locate the exact position of the two Dirac points in the ω-k space, the dispersion relation of the 3), which yields the following equation Here the wave vector k p = ω p /c.Note that the SPP dispersion relation requires that in the multilayer stack requires that the wave vector k x and k y must have real values.Hence, the only solution of Eq. ( 4) reads which stands for the location of one Dirac point in the ω-k space.Since Eq. ( 3) is an even function of the wave vector k y , there are two Dirac points symmetrically distributed along k y -direction with respect to the origin.It is interesting that Eq. ( 5) reveals that the positions of the Dirac points in the band structure are only related on the filling ratios f 1 and f 2 of the materials in the multilayer stack, rather than the actual thicknesses of different layers, which means that the positions of the Dirac points are determined once the filling ratios of the materials are specified.Moreover, the frequency with the frequency at which the effective permittivity ε (0) y = 0, as calculated from Eq. ( 2) based on the EMT.
It is necessary to investigate more about the dispersion relation around the Dirac points with respect to the wave vector k x and k y .Figure 2(a) shows the dispersion relation between the wave vector k x /k p and the frequency ω/ω p calculated from Eq. (3) with the wave vector with the constant ∆, a function of parameters ε ∞ , ε 2 , d 1 , and d 2 , as listed in the Appendix.
Here the constant ∆ has a value of 58.705.The linear dispersion relation between k x /k p and ω/ω p forms the cross section of the Dirac cone at the Dirac point in the ω-k x plane, which is similar to the case in the electron band structure of graphene, implying that the propagation of electromagnetic wave in the multilayer stack will be the optical analogue to the charge transport in graphene.
In addition, Fig.
under the condition of the wave vector k x /k p = 0, which implies two different bands as and because ε . Clearly, the two bands of the dispersion curves from Eq. ( 3) for the multilayer stack intersect at the Dirac point, and converge to the SPP dispersion curve when the wave vector k y /k p is increased, due to the surface plasmon resonance (SPR) in the condition of ε 1 = −ε 2 .On the contrary, the two bands of the dispersion curves based on the EMT of Eqs. ( 8) and ( 9) do not predict the SPR at large wave vector k y /k p , since the EMT does not take into account the layered structures of the multilayer stack, but they still intersect at the Dirac point.Furthermore, the dispersion relation from Eq. ( 3) for the multilayer stack also obeys a linear relation at the Dirac point, which is indicated by the black straight lines that obey the following equation .
However, the nonlocal ENZ point for the multilayer stack due to the optical nonlocality is shifted into the position of k x /k p = k y /k p = 0 but with a lower frequency, which is determined from Eq. (3).

III. ISO-FREQUENCY CONTOURS AND PROPAGATING MODES
The IFCs of multilayer stack represents the spatial dispersion in k-space, which are directly related to the propagating properties of electromagnetic waves.Displayed in Fig. 3 3) for the multilayer stack are plotted as red curves, while the IFCs obtained from the EMT of Eq. ( 7) are plotted as blue curves for comparison.It is found that the IFCs obtained from Eq. ( 3) consist of two branches at all frequencies, a hyperbolic-like branch and an ellipticlike branch, which is coincident with the band structure shown in Fig. 1(c), and the two branches join together at the Dirac point shown in Fig. 3(c).On the contrary, the EMT calculated IFCs consist of only one single branch, varying from hyperbolic-like to elliptic-like as the frequency is increased, and they possess similar curvatures as that of the IFCs from Eq. (3) when |k x /k p | 1, which is coincident with the long wavelength limit for the EMT.
It is shown in Fig. 3(c) that the EMT calculated IFC reduces into a straight line along the k y -axis at the frequency of the Dirac point, which is corresponding to the band structure described by Eq. ( 9).While the nonlocal IFC at the same frequency obtained from Eq. ( 3) shows a dramatic difference, revealing a giant optical nonlocality at the Dirac point.With the optical nonlocality, the permittivity will be a function of both the frequency and the wave vector, which can be analytically described in an approximate way based on a modification about the EMT of Eqs. ( 1) and ( 2), through the Taylor expansion of Eq. ( 3) under the conditions of x d 1 1, and k x d 2 1, which are all satisfied in the neighborhood of the Dirac point.By expanding Eq. ( 3) up to the third non-vanishing term, the approximated dispersion relation can be written in the form of Eq. ( 7) as with the modified nonlocal effective permittivity components of ε eff x = ε (0) x /(1 − δ x ), and , where the nonlocal modification factors δ x and δ y are the functions of both frequency and wave vector as and in which all the terms are listed in the Appendix.The approximated IFC based on Eq. ( 11) at the frequency of the Dirac point is plotted as black curves in Fig. 3(c), and it matches well with the IFC directly obtained from Eq. ( 3) near the Dirac point.The nonlocal modification factors clearly reveal that the modified nonlocal effective permittivity ε eff x and ε eff y are not only related with the frequency, but also strongly dependent on the wave vector components k x and k y .It is shown that at the Dirac point, the nonlocal modification factor δ x = 0, which x and the nonlocal effect does not affect the effective permittivity component vertical to the multilayer interface.However, there is significant difference for the modified nonlocal effective permittivity component ε eff y , where the nonlocal modification factor δ y = 1 at the Dirac point, leading to an indeterminate form of the by the giant optical nonlocality is also demonstrated.Finally, it is noted that although the study is carried out under lossless condition, the giant optical nonlocality still affects the propagation of the electromagnetic wave with a moderate material loss.

APPENDIX
The constant ∆ in the linear dispersion relation Eq. ( 6) between the wave vector k x /k p and the frequency ω/ω p can be calculated as follows in which .
The constants C 1 , C 2 , and C 3 in the linear dispersion relation Eq. (10) between the wave vector k y /k p and the frequency ω/ω p can be calculated as follows and with all terms listed as follows around the Dirac point.The band structure with the band 1 and the band 2 connected with the Dirac point is clearly illustrated.Furthermore, the linear dispersion at the Dirac point, a necessary condition for a Dirac point, is clear revealed by the Taylor expansion, and indicated by the black straight lines that obey the following equation 2(b) displays the dispersion relation between the wave vector k y /k p and the frequency ω/ω p based on Eq. (3) with the wave vector k x /k p = 0 around the Dirac point.For comparison, the SPP dispersion relation k y /k p = ε 1 ε 2 /(ε 1 + ε 2 )ω/ω p , and the dispersion relation based on the EMT are also plotted.It is noted that the EMT dispersion relation is calculated from the following equation

) The constants C 1 ,
C 2 , and C 3 are the functions of parameters ε ∞ , ε 2 , d 1 , and d 2 , as listed in the Appendix.Here the three constants have the values of C 1 = 58.705,C 2 = −20.0954,and C 3 = 0.282375.Equation (10) reveals the asymmetry of the dispersion relation between the wave vector k y /k p and the frequency ω/ω p with respect to the Dirac point, which leads to a titled cross section of the Dirac cone in the ω-k y plane.Additionally, it is noted that the location of the Dirac point is not at the center of the BZ, where the wave vector k x /k p = k y /k p = 0, due to the different periods of the multilayer stack along x-and y-directions.Therefore, the Dirac point cannot be mapped into an ENZ point in the ω-k space, although the frequency of the Dirac point is the same as the ENZ frequency predicted by the EMT of Eq. (2), where the effective permittivity ε (0) y = 0. Depicted in Fig.2(c), the EMT predicted ENZ point is located at the position of k x /k p = k y /k p = 0 and ω , the IFCs of five different frequencies near the Dirac point are plotted, including 636.577 THz [the nonlocal ENZ frequency from Eq. (3)], 641 THz (slightly below the Dirac point), 647.027THz [at the Dirac point from Eq. (5)], 651 THz (slightly above the Dirac point), and 671 THz (far above the Dirac point).The IFCs calculated from the band structure of Eq. ( ) at the Dirac point, which is coincident with the previous analysis that the Dirac point in the band structure of the multilayer stack cannot be mapped into the ENZ point.According to the IFCs with nonlocal effects calculated from Eq. (3), the corresponding eigenmodes and propagating patterns are analyzed for the electromagnetic wave propagating along y-direction in the multilayer stack in Fig.4, at three frequencies around the Dirac point including 641 THz, 647.027THz, and 651 THz, displayed in Fig. 4. The eigenmodes supported in one period of Al 2 O 3 -Au-Al 2 O 3 unit cell of the multilayer stack are illustrated with the distributions of the magnetic field amplitude H z and the absolute value of the magnetic field abs(H z ).When the frequency is below or above the frequency of the Dirac point, the propagating electromagnetic wave possesses two different eigenmodes, the symmetric mode and the anti-symmetric mode, with different propagating constants k y as marked on the IFCs in Fig. 3.The two eigenmodes degenerate into a single symmetric mode at the frequency of the Dirac point.Besides, an interchange of the two eigenmodes occurs when the frequency pass across the frequency of the Dirac point, due to a band inversion at the Dirac point as shown in Fig. 2(b).In addition, the giant optical nonlocality near the Dirac point of the multilayer stack can result in a unique optical phenomenon in the propagation of the electromagnetic wave.The propagation of a TM polarized Gaussian beam (with E x , H z , and k y ) is considered at the three different frequencies, and the distributions of the absolute value of the magnetic field abs(H z ) are shown in Fig.4.For comparison, both the multilayer stack and the corresponding effective medium are simulated.Regarding the normal incidence, it is found that the Gaussian beam is scattered into the similar diverging patterns in both multilayer stack and effective medium, when the frequency is below and above the frequency of the Dirac point, due to the common sharp curvature of the IFCs near the Dirac point in Figs.3(b) and 3(d).However, at the frequency of the Dirac point, the propagation of the electromagnetic wave is extraordinary.In the multilayer stack, due to the degeneration of the symmetric mode and the anti-symmetric mode, the joint of the two IFC branches at the Dirac point flatten the sharp curvature, leading to a splitting of the Gaussian beam into two mirrored propagating directions, as shown in Fig.4(b).The beam splitting phenomenon in the multilayer stack due to the optical nonlocality can be looked as the optical analogue to the giant nonlocality enhanced long-range flavor current splitting in the quantum Hall effect of graphene at the Dirac point.On the contrary, the incident electromagnetic wave is totally prevented from propagating into the corresponding effective medium due to the large impedance mismatch, which is coincident with the EMT calculated IFC shown in Fig.3(c).IV.CONCLUSIONSThe giant optical nonlocality near the Dirac point for the multilayer stack have been revealed in lossless metal-dielectric multilayer metamaterials by studying the band structure in the three dimensional ω-k space based on the transfer-matrix method with the optical nonlocal effect.The exact location of the Dirac point in the band structure of the multilayer stack is determined and the dispersion relation around the Dirac point is investigated in details.Based on the mode analysis, it is proved that the degeneracy of the symmetric mode and the asymmetric mode of the coupled SPPs is the origin of the Dirac point.Meanwhile, the position shift of the ENZ point in the band structure affected by the giant optical nonlocality near the Dirac point is also explored, and a nonlocal modification on the dispersion relation based on the EMT including the optical nonlocality near the Dirac point is derived.Furthermore, the giant optical nonlocality near the Dirac point is revealed by means of the IFC analysis, and the extraordinary beam splitting at the Dirac point induced

FIG. 2 .FIG. 3 .FIG. 4 .FIG. 1 .FIG. 2 .FIG. 3 .
FIG. 2. (a) The dispersion relation k x /k p ∼ ω/ω p based on Eq. (3), under the condition of k y/k p = ε 2 d 1 d 2 / ((d 2 − d 1 )(ε ∞ d 1 + ε 2 d 2 ))in red curves.The band 1 and band 2 are connected at the Dirac point, located at the position of k x /k p = 0 according to Eq. (5).The black lines indicate the linear dispersion in the neighborhood of the Dirac point consistent with Eq. (6).(b) The dispersion relation k y /k p ∼ ω/ω p near the Dirac point based on Eq. (3), with respect to k x /k p = 0 in red curves.The dispersion relation obtained from the EMT and the SPP dispersion relation are plotted in dot-dashed black curve and dashed blue curves, respectively.All the dispersion curves intersect at the Dirac point.The dispersion curves obtained from the multilayer stack converge to the SPP dispersion when the wave vector k y /k p increases, due to the SPR.The linear dispersion relation in the neighborhood of the Dirac point is plotted in black lines based on Eq. (10).(c) The positions of the ENZ determined from the EMT and the multilayer stack including the optical nonlocality are marked for comparison.