Quasi-confined ENZ mode in an anisotropic uniaxial thin slab.

In this paper, we present a numerical modal study of a simple slab, made of an uniaxial anisotropic material having an "epsilon-near-zero" (ENZ) dielectric function, surrounded by vacuum. We use two Drude models with a different plasma frequency for the direction parallel and perpendicular to the slab surface as toy models to study the effect of uniaxial anisotropy of type I (∊‖ > 0, ∊⊥ < 0) and type II (∊‖ < 0, ∊⊥ > 0) on the different electromagnetic modes of the system. In addition to the so-called ENZ mode, studied in detail by Campione et. al [ Phys. Rev. B91, 121408(R) (2015)], the slab can support quasi-confined (QC) mode in the type I and type II anisotropy frequency ranges. We show that those modes exhibit a strong electric field enhancement, caused by the ENZ character of the dielectric function. In strong contrast with the ENZ mode, QC modes can have a strong electric field enhancement for thick slabs, with a Fabry-Perot-like electromagnetic field distribution spanning over the whole slab thickness. This opens the way for large electric field enhancement in thick slabs with QC ENZ modes. Thick slabs also allow metamaterial designs, giving the possibility to engineer the anisotropy of the effective dielectric function, opening interesting perspectives for the control of field enhancement of the ENZ QC modes and their integration in operating devices, such as detectors, sources, or modulators.


Introduction
In their pioneering work [1], Engheta and co-workers predicted the tunneling of electromagnetic energy through narrow channels with the use of material having a dielectric function close to zero (epsilon-near-zero, or ENZ). Experimental demonstrations came shortly after, in the microwave regime [2][3][4]. These works triggered investigations of the peculiar properties of materials with ENZ behavior. Many different materials and many different effect have been reported, either theoretically or experimentally, such as polarization control [5], control of the radiation phase pattern of emitter close to ENZ structures [6], boosting of non-linearities [7,8], light modulators [9], perfect absorber [10][11][12][13], and thermal emitter [14] to name a few. We refer the reader to a recent review for more informations on this topic [15]. ENZ materials can be found naturally in polar materials close to optical phonon resonances, in transparent conductive oxide close to their plasma frequency, but have also been engineered to reach visible wavelength using metamaterials [16] or organic molecules [17]. Strategies to enhance the ENZ properties by loss-compensation have also been proposed [18].
Another important topic is the electromagnetic modes supported by thin slabs of ENZ materials. The so-called ENZ mode was shown to be the evanescent counterpart of the Berreman mode [19] and was used as a first application to build a THz reflectivity modulator [20]. The mode is characterized by a strong enhancement and confinement of the electric field perpendicular to the slab surfaces. Strong coupling with other type of resonances has been experimentally observed [21][22][23][24]. A detailed theoretical study of the ENZ mode was recently published [25] where it was shown that this ENZ mode shows field enhancement only for a limited range of slab thickness and wavevector. In particular a rule of thumb was given for the slab thickness limit: d λ p /50, where λ p is the wavelength corresponding to the plasma frequency for a simple Drude model. This appears as a significant limitation when using the mode to enhance the electrical field for absorption or non linear effect, as the volume of the active material is significantly reduced.
Here, we extend this theoretical work by considering an uniaxial anisotropic slab with ENZ properties. We focus on type I ( > 0, ⊥ < 0) and type II ( < 0, ⊥ > 0) anisotropy. Both types give birth to quasi-confined (QC) modes whose dispersion and spatial structure strongly differs from that of the ENZ mode. We show in this paper that they have a large field enhancement due to the ENZ properties of the slab. In strong contrast with the ENZ mode, the field enhancement can be maintained even at large slab thicknesses, and the electric field spatial distribution spans over the whole slab thickness.
The paper is organized as follows: in section 2 we define the geometry and model describing the system under study, in section 3 we discuss the QC mode and their field enhancement for different slab thicknesses, in section 4, we discuss the coupling of these modes. Finally in section 5 we briefly discuss possible experimental realization of such anisotropic ENZ materials.

Theory
The system under study is a slab of thickness d made of an anisotropic uniaxial material. The geometry is illustrated in Fig. 1. The surfaces are perpendicular to the z direction and lie in the (x, y) plane. Both interfaces are respectively placed at z = 0 and z = d. The anisotropic material is described by a uniaxial dielectric tensor: (ω).
The subscript ⊥ and are taken relatively to the (x, y) plane: ⊥ describes the dielectric function along the z direction ( Fig. 1). We introduce a simple Drude model as toy model for the slab permittivity in both directions and consider vacuum as the dielectric surrounding the slab: We note that for this model, when ω = ω p 1 − γ 2 ω 2 p ≈ ω p , the real part of the dielectric function is zero. As explained in [20] and [25] the electric field along the perpendicular direction (E z ) is then strongly enhanced. This comes from the continuity of the electric displacement at the interface between the slab and vacuum. By using constitutive equations we can express the electric field along the z direction inside the slab (E zs ) as a function of the dielectric function of the slab ( ⊥s ) and the electric field outside the slab (E zv ): We introduced the K E N Z factor that represents the intensity enhancement inside the slab due to the ENZ behavior of the dielectric function. From Eq. (4) one can clearly check that if the dielectric function of the slab ⊥s goes to zero, E zs is diverges. In reality, because of losses, ⊥s is never strictly zero, and the enhancement, given by |1/ ( ⊥,s )| 2 , is finite. The lower the imaginary part is, the stronger the enhancement is. We will call the frequency range where K E N Z > 20 the ENZ range. This value corresponds to ω p ± 10% as discussed in [25]. The dielectric function ⊥s , K E N Z factor and the corresponding ENZ range are represented in Fig. 2 for a Drude model with parameters following [25]: ω p⊥ =10000 cm −1 , γ ⊥ = 100 cm −1 . Fig. 2. ENZ frequency range and enhancement associated with a Drude model describing the dielectric properties of a slab in the direction perpendicular to the interfaces ( ⊥s ). Left: real and imaginary part for a dielectric function described by a Drude model with ω p⊥ =10000 cm −1 , γ ⊥ =100 cm −1 . Right: corresponding K E N Z factor, representing the enhancement of the electric field intensity. In both figures, the green shaded area corresponds to the ENZ range defined as K E N Z > 20.
However, as pointed out in [25], the ENZ enhancement is limited to thin slabs. In order to overcome this limitation, we now move to anisotropic materials. We study transverse plane waves propagating along the x direction and polarized with the electric field along the (x, z) plane (TM polarization). Using Maxwell's equations we derive the dispersion relation for the anisotropic slab (see the Appendix for derivation, with the implicit exp (−iωt) dependence): where k ⊥v and k ⊥s are the wavevector component along the z axis respectively in vacuum and in the slab, given by Eqs. (6) and (7): with a square root determination given by (k ⊥v ) + (k ⊥v ) > 0. The anisotropy is introduced by setting ω p to a different value than ω p⊥ . Between roughly ω p and ω p⊥ , the real part of s and ⊥s will have opposite signs. This is illustrated in the left panel of Figs. 3 and 4. In this situation, with ⊥s close to zero, we can see from Eq. (7) that k ⊥s will be real and will take large values. For any value of the slab thickness d, one can always find an integer p such that 2k ⊥s d = p2π, because k ⊥s can take any large value as ⊥s tends to zero. Hence, the electromagnetic modes in this frequency range will have Fabry-Perot-like spatial field distribution along the z-direction inside the slab. These modes are the QC modes that we will study in the next section. We will show that using such modes, strong electric field enhancement can be achieved even for large slab thicknesses.
Finally, throughout this paper and in the Appendix, we will refer to symmetric and antisymmetric modes considering the electric field perpendicular to the slab (E z ), as historically used by Sarid [26] in his seminal paper on long-range (LR) and short range (SR) mode. With this convention, the LR mode is symmetric, while the SR mode is anti-symmetric. The anisotropy changes slightly their dispersion relation (see the Appendix for more details).

Quasi-confined ENZ mode
The QC modes were mostly discussed in the framework of electron-phonon interaction in GaN/AlGaN quantum heterostructure (see [27] and citing articles) and more recently in the context of hyperbolic materials [28,29] and quantum transport in graphene supported by h-BN [30] where they are referred to as Hyperbolic Phonon Polaritons (HPP). In both cases, the anisotropy is naturally present thanks to the optical phonons and wurzite structure of GaN and h-BN. However, these modes are always calculated without taking losses into account. Furthermore, the electric field enhancement due to the ENZ properties of the materials was not discussed so far.
We plot in the left panels of Figs. 3 and 4 the dielectric function used in our system respectively for type I (ω p < ω p⊥ ) and type II (ω p > ω p⊥ ) anisotropy. The ENZ range (where K E N Z > 20) is represented by a green shaded area, while the frequency range where QC modes exist is represented by a hatched area.
The solutions for the modal Eq. (5) are numerically found in the complex plane, using two different algorithms (a home made Newton algorithm and Reticolo [31]), both yielding similar results. We choose here a representation using complex frequency and real wavevector scheme, as it is relevant to Local Density of Optical States (LDOS) and pulsed excitation experiments, while real frequency and complex wavevector is more suited to propagation of fields. See [32] for a detailed discussion We plot in the right panels of Figs. 3 and 4 the corresponding dispersion relation of the modes supported by a slab of 5 nm thickness. For clarity, only QC1 and QC2 are shown. In reality there is theoretically an infinity of such modes of higher order (see the Appendix, Fig. 12 for the type I case).
The symmetric LR (black continuous line) enters the ENZ range for a limited wave-vector range, where the mode gets then a strong E z enhancement, as in [25]. The QC modes have ω p as their asymptotic frequency at large k x , and they tend interestingly towards ω p⊥ at small wave-vector. The first order QC mode, symmetric and anti-symmetric are the most dispersive, higher order disperse less and less as the mode number increases.
We can clearly see that for this level of anisotropy and this slab thickness the QC modes mostly stay in the ENZ range. The spatial electromagnetic field distribution of the different modes are shown in Fig. 5 for the case of type II anisotropy. Results for the type I are very similar and are not shown here. As we get complex frequency from our modal search, we can extract the quality factor of the modes (see Fig. 16 in the Appendix). QC modes have a slightly higher quality factor than the LR and SR modes, but are of the same order of magnitude. For the case described in Fig. 4, a quality factor of the order of 100 is found for QC modes, corresponding to a lifetime of roughly 100 fs. The field distribution of the QC modes strongly differs from the LR and SR modes. The spatial structure of |E z | for the QC1 shows two nodes, and 4 for QC2. QC modes of higher order will show more nodes. Note however that the field of higher modes varies over distances so small that the validity of the permittivity model may fail for real materials. The fields outside the slab have an exponential decay, typical of surface mode, and very similar to the LR and SR modes. As can be seen in the upper right panel of Fig. 5, both QC modes have a larger enhancement than the SR mode: the maximum of |E z | is 73 times larger for QC1 and 296 times larger for QC1 than for the LR ENZ mode. This is simply due to the fact that the QC modes are closer to the frequency where K E N Z is maximum (≈ ω p⊥ =10000 cm −1 ) than the LR ENZ mode. Note that both symmetric and anti-symmetric QC mode are in the ENZ range and benefit from the strong E z enhancement.
From [25], we know that if we increase the slab thickness further than λ p /50 ( 20nm in our case), the symmetric LR mode dispersion relation does not reach the ENZ range and thus loses its E z field enhancement. For the QC modes, things are different. The dispersion of the QC mode is stronger as the thickness of the slab gets larger, but they still tend to ω p⊥ at small wavevector. Furthermore, the QC range is located roughly between ω p⊥ and ω p , and the ENZ range is centered around ω p⊥ . If the anisotropy is weak (ω p close to ω p⊥ ), then the QC range can be completely inside the ENZ range. This is a way to control the dispersion of the QC mode and to keep them in the ENZ range. This can be seen in Fig. 6 where we show the dispersion relation of the modes supported by a 100nm thick slab with ω p = 1.05ω p⊥ . The LR and SR modes do not enter the ENZ range, while the QC mode are all located in the ENZ range. Fig. 6. Dispersion relations for a slab of 100nm thickness. Solid lines indicate symmetric modes and dashed lines anti-symmetric modes. The green shaded area correspond to the ENZ range, and the hatched area correspond to the QC mode range. The LR mode is not ENZ anymore. ω p⊥ = 10000cm −1 , ω p = 1.05ω p⊥ , γ = γ ⊥ = 100 cm −1 , d = 100 nm.
The fields distribution of different modes are shown in Fig. 7, for k x = 5 × 10 7 m −1 .The LR mode does not have any enhancement of |E z | (top right in Fig. 7) as it is not in the ENZ range. The electromagnetic field profile of the LR and SR are the "usual" distributions with exponential decays inside and outside the slab. On the other hand, the QC mode are located entirely in the ENZ range. They show a clear |E z | enhancement and their spatial field distribution spans along the whole thickness of the slab with the same structure as in the thinner slab. Note that in this configuration, QC1 is the furthest away from the ENZ frequency (ω p⊥ ), meaning all higher order QC modes have a larger enhancement. This is seen for QC2 in Fig. 7, but higher order QC modes are not represented for clarity reasons (see Appendix, Fig. 13). Higher losses change only very slightly the dispersion relations.
To summarize our findings, QC mode can show a strong E z enhancement due to the ENZ properties of the slab material when their dispersion relation enters the ENZ range. This happens naturally at small wavevector and can be controlled by the level of anisotropy. Remarkably, a small anisotropy (ω p⊥ close to ω p in our case) keeps the QC mode in the ENZ range. The obtained enhancements are larger than the classical ENZ mode as the modes frequencies come closer to the frequency where K E N Z is maximum. The electromagnetic field distribution shows complex structures spanning over the whole slab thickness, even for very large slabs. More importantly, QC ENZ mode can still exist for thick slabs. This make them good candidates for the design of absorbers, or efficient non linear materials, as the interaction volume can be large while still benefiting of a strong field intensity.

Quasi-confined mode excitation
In this section, we briefly discuss how to excite the QC mode. First, we consider a dipolar point-like source located at 10 nm above a 200 nm slab of anisotropic materials (ω p⊥ = 10000cm −1 , ω p = 1.4ω p⊥ , γ = γ ⊥ = 100 cm −1 ), with different dipole orientations. Figure 8(a) shows the Purcell factor (P/P vac , where P is the power emitted by the dipole above the slab, and P vac the power emitted by the dipole in vacuum). We can see a strong enhancement around 8000 cm −1 corresponding to the LR mode of the slab. An enhancement of around one order of magnitude can be seen between ω p⊥ and ω p , corresponding to the QC modes. This limited enhancement, despite the large number of modes can be explained by the very strong confinement of the QC modes inside the slab, leading to a poor overlap of the dipole and QC modes field.
Then, we turn to the classical Kretschmann configuration, using a silicon prism (n Si = 3.5). Figure 8(b) shows the absorption as a function of frequency and wavevector. Between ω p⊥ and ω p ,we can clearly see the QC modes, that still exist in this non-symmetrical dielectric environment. The coupling in this configuration is strong but the absorption clearly decreases as the frequency goes to ω p⊥ = 10000cm −1 . The vacuum light line is represented as a green line. The QC modes still exist inside the light cone, which means that they could be directly excited by vacuum propagating light. This can be seen in Fig. 15 in the Appendix.
We represent in Fig. 8(c) the spatial distribution of |E z | as a function of z and frequency for the Kretschmann configuration, for an incident angle of 50 degrees. This angle corresponds to the dashed red line in Fig. 8(b). The structure of the QC modes in this configuration looks like the one from Figs. 5 and 7. The mode at 13000 cm −1 is the QC 1 anti-symmetric mode, having two maximums at the edge of the slab. The mode at 11500 cm −1 is the QC 1 symmetric mode, having two maximums at the edge of the slab and one in the middle, and so on. The corresponding |E x | field are illustrated in Fig. 14 in the Appendix. We note that QC modes have been observed by other groups with other methods. Signature of electrical excitation of QC modes have been reported for electrons travelling in a graphene sheet supported a h-BN thin film [30]. Direct observation of QC modes have been reported by apertureless scattering interferometric near-field optical microscopy on a 150 nm thick h-BN film [33].

Materials with anisotropy
In this last section we discuss briefly the materials that exhibit uniaxial anisotropy and could be used to design active devices like optical modulator, detectors or sources. In particular, the enhanced electric field in the direction perpendicular to the slab interface is of great interest for interactions with intersubband transitions in semiconductor quantum wells as they interact only with this component of the electric field. To have QC mode with strong electric field enhancement, we need a combination of anisotropy of type I or type II, and an ENZ dielectric function in the z direction.

Natural materials
Both of these conditions are naturally found in crystalline wurzite systems (GaN, AlGaN, h-BN, SiC-4H and SiC-6H, sapphire) close to their longitudinal optical phonon frequency. They usually lie in the far infrared (wavelength > 7µm). At the LO phonon frequency in the z direction, the dielectric function is ENZ with low losses (depending on the crystal quality). Because of the anisotropy in the crystalline structure, the LO phonon in the (x,y) direction is shifted compared to the z direction and can lead to type I or type II anisotropy.
As discussed in this paper, and experimentally shown in numerous publications, free electrons in doped semiconductors can lead to ENZ dielectric function. Anisotropy in the crystalline structure can lead to an anisotropy in the effective mass of the electrons, leading to different plasma frequencies in different directions.

Quantum materials
The dielectric function of a quantum well (QW) with intersubband transitions can be described by an anisotropic tensor, where the electrons in the well are described by a Drude model in the direction perpendicular to the QW growth direction, and by a resonant term in the growth direction describing intersubband transitions. These transitions are well studied, in particular in the GaAs/AlGaAs system. We also note that using strongly doped large QW, collective behavior of electrons can lead to a redistribution of the energy levels and generate a strong transition [34]. The condition to have ENZ behavior is to have a transition strong enough so that it can overcome the matrix static dielectric function to effectively obtain a negative real part of the dielectric function, and a long lifetimes (low losses), to obtain good ENZ properties.

Artificial metamaterials
The important conclusion of this study is that the ENZ-enhanced QC modes exist for large slab thicknesses, and have an electric field that spans along the whole slab. This implies that the mode can "see" a metamaterial as an effective medium. Lamellar metamaterials, also known as hyperbolic materials can be described with the help of an uniaxial anisotropic effective dielectric function. An anisotropy of type I or type II, as well as a designed ENZ wavelength, can be achieved by combination of metals and dielectrics, transparent conductive oxide and dielectrics, or doped and un-doped semiconductors.

Conclusion
In this paper we used a simple toy model to study the electromagnetic modes supported by a slab of anisotropic material showing type I or type II anisotropy and ENZ properties. Both types of anisotropy lead to the appearance of QC modes. We showed that these modes can present strong electric field enhancement due to the ENZ properties of the slab. Even for thick slabs, the spatial electromagnetic field distributions of the QC modes spans over the whole thickness of the slab, with strong electric field enhancement typical of ENZ modes. This features opens the way for metamaterial designs, where the level of anisotropy can be tuned to keep the whole QC modes dispersion relation in the ENZ range. An important remaining question is the coupling of the QC modes to propagative light. Our calculation show that some modes stop at the light line, while other cross it, with a small discontinuity in the imaginary part of the mode. Our understanding of this behavior is limited to date and will be the subject of further study, as well as other coupling schemes.

Dispersion relation derivation
We recall the Maxwell equation without sources and external charges: We study transverse plane waves propagating along the x direction and polarized with the electric field along the xz plane (TM polarization) We look for exponentially decaying fields away from the slab interfaces. Using convention from Fig. 1, we have: for the field in above the slab (z<0). Below the slab (z>d), we have: Inside the slab, the fields are described by two exponentially decaying fields from both interfaces: Using Ampere's law (Eq. (11)), we can express E x : This gives explicitly, at time t=0 and position x=0: From the continuity of H y and E x at z=0 and z=d we get the following system of equations: By inserting Eq. (21) into Eq. (22), we end up with the following relation:

Effect of anisotropy on LR and SR modes
The global behavior of the LR and SR mode is similar to the isotropic case, except for a change of the asymptotic frequency. This is shown in Fig. 9. Âń In particular for a strong anisotropy, the asymptotic frequency of the LR and SR modes enters the ENZ frequency range. The asymptotic frequency can be determined for the expression of k in the case of a single interface between air and the anisotropic material [35]: The asymptotic frequency correspond to the frequency when k , i.e. when 1 − ⊥ = 0. When we inject our Drude model in this expression we find that the asymptotic frequency ω asympt is given by: where losses were neglected for simplicity. From Eq. (25), it appears that ω asympt will never reach ω p⊥ (maximum ENZ). However,we can try to get ω asympt as close as possible to ω p⊥ . For this we can write ω asympt = aω p⊥ in Eq. (25), and we get and expression for ω p : The function clearly tends to infinity for a=1, but for a=0.9, one gets ω p = 2ω p⊥ , and for a=0.99, ω p = 7ω p⊥ . These values are high but sufficient for the asymptotic frequency to be in the ENZ range. This implies that the symmetric LR mode will be in the ENZ range, independently of its wavevector. This is in strong contrast with the isotropic case where the mode enters the ENZ range only for a limited wavevector range as described in [25]. If the slab thickness is increased, the dispersion relation of the symmetric and anti-symmetric modes reach their asymptotic value for smaller wavevector. It is then possible to bring even the anti-symmetric mode in the ENZ frequency range, where E z is strongly enhanced (see Fig. 10). However this comes at the cost of strongly confined electric fields, with decay length so short inside the slab that the Drude models used in our calculation might fail (see Fig. 11). The strong anisotropy requirements, and the electromagnetic field distribution confined at the slab surfaces with decays close to the limit of locality for the definition of the dielectric function makes the LR anisotropic ENZ mode not so attractive for applications. Fig. 11. Electromagnetic field distribution for the LR symmetric mode (top row) and SR asymmetric modes (bottom row) for different levels of type II anisotropy. ω p⊥ = 10000cm −1 , γ = γ ⊥ = 100 cm −1 , d = 2 nm.

Dispersion relation in color plot for a thin slab
See the colormap plotted in Fig. 12. The quantity plotted is: where DR is the left part of Eq. (5), with ± = + for symmetric modes (left) and ± = − for anti-symmetric modes (rigth). Losses have been reduced to 10 cm −1 for clarity. Modes appear as yellow lines.

Dispersion relation in color plot for a thick slab
See Fig. 13. Same comments as the previous section.

Coupling to QC modes
In Fig. 14 we represent the spatial distribution along the z direction of both |E x | and |E z | as a function of frequency, at an angle of 50 degree, for a 200nm anisotropic slab on a silicon prism, with ω p⊥ = 10000cm −1 , ω p = 1.4ω p⊥ γ = γ ⊥ = 10 cm −1 , d = 200 nm. This corresponds to Figs. 8(b) and 8(c) in the article. In Fig. 15 we present the absorption for an anisotropic slab surrounded by vacuum, as well as the spatial distribution along the z direction of both |E x | and |E z | as a function of frequency, at an angle of 50 degree, with ω p⊥ = 10000cm −1 , ω p = 1.4ω p⊥ γ = γ ⊥ = 10 cm −1 , d = 200 nm.