Waveguide Modes in Weyl Semimetals with Tilted Dirac Cones

We theoretically study unattenuated electromagnetic guided wave modes in centrosymmetric Weyl semimetal layered systems. By solving Maxwell's equations for the electromagnetic fields and using the appropriate boundary conditions, we derive dispersion relations for propagating modes in a finite-sized Weyl semimetal. Our findings reveal that for ultrathin structures, and proper Weyl cones tilts, extremely localized guided waves can propagate along the semimetal interface over a certain range of frequencies. This follows from the anisotropic nature of the semimetal where the diagonal components of the permittivity can exhibit a tunable epsilon-near-zero response. From the dispersion diagrams, we determine experimentally accessible regimes that lead to high energy-density confinement in the Weyl semimetal layer. Furthermore, we show that the net system power can vanish all together, depending on the Weyl cone tilt and frequency of the electromagnetic wave.These effects are seen in the energy transport velocity, which demonstrates a substantial slowdown in the propagation of electromagnetic energy near critical points of the dispersion diagrams. Our results can provide guidelines in designing Weyl semimetal waveguides that can offer efficient control in the velocity and direction of energy flow.


I. INTRODUCTION
With the recent discoveries of a number of Weyl 1 semimetal (WS) compounds 2- 18 and their intriguing optical properties 4,[19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36] , there has been a surge in research that investigates their potential roles in future optical systems. When studying the electromagnetic (EM) properties of WS systems, it is important to take into consideration their inherent anisotropic nature. Depending on the polarization of the EM wave relative to the crystalline axes, a WS can be tuned to exhibit a metallic or insulating type of EM response, as described by its permittivity tensor . This is reflected in the real parts of the appropriate components of changing sign, or vanishing altogether, as occurs in epsilon-near-zero (ENZ) materials. [37][38][39][40][41][42][43][44][45][46][47][48][49] Various peculiar effects can arise in materials exhibiting an ENZ response 46 , including light tunneling, vanishing group velocity, and perfect absorption. 43 By exploiting the anisotropic tunable EM response of Weyl semimetals, new classes of optical devices, such as waveguides, can be designed that afford dissipationless, highly localized energy transport with minimal attenuation and signal distortion.
For anisotropic waveguides, the relative signs of the components of can also be crucial for the optimization of energy transport. For example, a uniaxial anisotropic system that has permittivities of opposite sign can have an isofrequency surface in the form of an open hyperboloid. Such hyperbolic metamaterials have been discussed in the context of waveguides 52 to generate forward and backward surface modes, 50 and to enhance EM fields. 51 Thus, the anisotropic EM response of a WS must be accounted for when discussing its role as a waveguide in the ENZ regime. Indeed, the sign differences that can arise between the diagonal permittivity components, can afford greater control and increased subwavelength EM field confinement.
In addition to the dielectric response of Weyl semimetals, the tilting of the Weyl cones can induce observable effects in the conductivity and polarization. 32 By reversing the conical tilt direction, the right-hand and left-hand responses of the Weyl semimetal become reversed, which can affect the absorption of circularly polarized light. 29 It was found that the chirality or sign of the tilt parameter in the absence of time-reversal and inversion symmetries can change the sign of the Weyl contributions to the absorptive Hall conductivity 33 . For Weyl semimetal structures much thinner than the incident wavelength, EM radiation can be perfectly absorbed for nearly any incident angle, depending on the tilt of the Weyl cones and chemical potential. 48 Moreover, the broad tunability of the chemical potential in a WS makes it a promising material for photonics and plasmonics applications 18,20,[24][25][26][54][55][56][57][58] . The chiral anomaly in a WS can alter surface plasmons and the EM response [24][25][26][55][56][57][58] . It has been shown theoretically 26 that measurements of the optical conductivity and the temperature dependence of the free carrier response in pyrochlore Eu 2 Ir 2 O 7 is consistent with the WS phase, and that the interband optical conductivity reduces to zero in a continuous fashion at low frequencies as predicted for a WS. The surface magnetoplasmons of a Weyl semimetal can turn to low-loss localized guided modes when two crystals of the WSs with different magnetization orientations are connected. 24 In this paper, we investigate guided wave modes that can arise in an anisotropic Weyl semimetal waveguide structure. We consider a finite-sized Weyl semimetal layer on top of a perfectly conducting substrate. Our study includes both analytic and numerical results that determine the EM fields and characteristics of the energy flow in WS waveguide systems. We consider broad ranges of conical tilting, and EM wave frequencies. We show that by appropriately tuning the frequency and conical tilt, both of the diagonal components of can achieve a lossless ENZ response. From the obtained dispersion relations, which reveal the permitted EM solutions, we construct the EM fields and identify regions of high energy density confinement in the WS layer. We demonstrate controllable net system power that can be made to vanish at critical points of the dispersion curves, where the propagation of electromagnetic energy comes to a standstill. Our findings offer systematic ways to design waveguide structures that can continuously reduce the energy velocity from a fraction of the speed of light to effectively zero.
The paper is organized as follows: In Sec. II A, we summarize the dielectric response of the WS by presenting each component of the permittivity tensor, and give details of the approximations made. Next, in Sec. II B, we solve Maxwell's equations and obtain the dispersion relations for the corresponding guided electromagnetic wave modes. We also define the averaged energy density, Poynting vector, and velocity of energy flow. In Sec. II C, we examine the obtained analytical expressions numerically and discuss our main results. Finally in Sec. III, we give concluding remarks.

II. MODEL AND RESULTS
In this section, we outline the salient features of the electromagnetic response of a WS by providing analytical expressions for the components of the permittivity tensor, as well as the numerical approach used. The Weyl semimetal is modeled as a system with broken time-reversal symmetry and two Weyl nodes. This model can be achieved through the stacking of multiple thin films involving a topological insulator and ferromagnet blocks. 4

A. Permittivity and approximations
The tilt of the Weyl cones are described by the parameters β ± , and for centrosymmetric materials with broken time reversal symmetry, we apply the condition β = β + = −β − . The separation between two Weyl points in the z direction in momentum space is defined by 2|Q|, where the sign of Q depends on the sign of the magnetization.
The permittivity tensor for the WS takes the following gyrotropic form: 48 where (ω) represents the equal components of (ω) parallel to the interfaces, i.e., xx (ω) = yy (ω) = (ω). Here the components are normalized by 0 , the permittivity of free space. The (ω) components can be written analytically as, 48 where µ is the chemical potential, with µ ≥ 0, α = e 2 /(4π 0 v F ), and Θ(X) represents the usual step function. Generally, the energy cut-off Γ is a function of the tilt parameter. Nevertheless, in our calculations, we choose a large enough cut-off and neglect the contribution of β to Γ. It is evident that (ω) are independent of the tilting parameters β ± . For the situation when the Weyl cones are not tilted (β ± = 0), we have zz (ω) = (ω), and the offdiagonal frequency dependent component γ(ω) reduces to, γ(ω) = 2v F αQ(πω) −1 . The imaginary term in Eq. (2) describes the interband contribution to the optical conductivity, which exists only when the frequency ω of the EM wave satisfies ω > 2µ. For the guided waves of interest, we are interested in the situation where there is no dissipation in the medium. As Eq. (2) shows, the range of frequencies in which is real and positive corresponds to Similarly, the permittivity component zz depends on frequency and chemical potential, however it also depends on the tilt β ± . The zz component involves a complicated integral, and its derivation is presented elsewhere. 48 For frequencies satisfying there is no dissipation (for all diagonal components), and zz can be written analytically as, 48 which is valid for |β ± | < 1. Below, we show that over the frequencies of consideration, zz and can be tuned to achieve an ENZ response at different ω. We also show that there are a range of frequencies where and zz can have opposite signs. This signs of the orthogonal permittivity components will be shown below to strongly affect the types of modes that will propagate along the waveguide interface.
The remaining off-diagonal term γ can play an important role in changing the polarization state of electromagnetic waves interacting with the WS via Faraday and Kerr rotations. 23,27 Variations in the gyrotropic term can also cause shifts in the surface plasmon frequency. 59 It can be written 1. (Color online). Schematic of the waveguide configuration involving a Weyl semimetal (in region 1) with width d on top of a perfectly conducting substrate. In this configuration, the wavevector k resides in the plane of the WS (along x), while the separation of Weyl nodes is perpendicular to the layers (along z). The surrounding medium (region 0) is taken to be vacuum. analytically in the limit ω → 0: Thus, for fixed β, and small tilting, γ is a linear function of µ, declining as the chemical potential increases. For arbitrary β and ω, it is necessary to resort to numerics to determine γ (see Appendix).
We now make use of the expressions above governing the relevant components of to demonstrate that by properly tailoring the EM response, the waveguide structure can effectively transport localized EM energy over a broad range of frequencies and tilting of the Weyl cones ( 0 < β < 1). When characterizing the nontrivial behavior of in the WS, there are several relevant parameters to consider, including the chemical potential, frequency of the EM wave, tilt of the Weyl cones, and the node separation parameter Q (which is taken to be positive). For simplicity, we consider here the configuration where β + = −β − = β.

B. Guided wave modes and energy confinement
We now investigate the EM modes for the configuration shown in Fig. 1, which consists of a planar Weyl semimetal (region 1) adjacent to a metallic substrate with perfect conductivity (PEC). The propagation constant k x is invariant across each layer. For both regions 0 and 1, we implement Maxwell's equations for time harmonic fields, where i = 0 or 1. Within the WS, the propagation vector k 1 replaces the spatial derivatives, transforming Maxwell's equations into the forms, k 1 × E 1 = ωµ 0 H 1 and k 1 × H 1 = −ω 0 E 1 . These two equations together result in the following expression for the E 1 field in k-space: Using k 1 = k xx + k 1zẑ , and the identity where k 0 = ω/c, and k ⊥ = k 2 1z + k 2 x . Due to translational invariance, we also have k 1x = k x . The presence of all three E field components in Eq. (9) is a consequence of the centrosymmetrc nature of the WS that can couple different polarization states, depending on the WS material parameters and layer width.
Taking the determinant of the matrix in Eq. (9) and setting it equal to zero, gives the dispersion equation for the wavevector k 1z in the semimetal: Solving Eq. (10) for k 1z results in two types of solutions denoted by k + and k − : The dispersion equation (10) can now be compactly written in terms of the two types of waves: When considering the transmission of EM surface waves, the traveling wave is confined to be within, or adjacent to, the waveguide walls. For the configuration shown in Fig. 1, where the x -y plane is translationally invariant, we consequently search for modes that are localized around the interface separating the WS and vacuum. As with all open waveguide structures, the WS waveguide spectrum consists of a finite discrete set of guided modes with purely real longitudinal propagation constants and an infinite continuum of radiation modes. Although outside the scope of this work, leaky waves can also arise, which have discrete complex longitudinal propagation constant solutions to the dispersion equation, and that decay longitudinally but do not obey the transverse radiation condition.
The magnetic field components in the vacuum region, H 0 , are thus written accordingly: where k x is the propagation constant, k 0z = k 2 x − k 2 0 , and the coefficients r 1,2,3 are to be determined upon application of the boundary conditions. Using the Maxwell equation ∇·H 0 = 0, a simple relation is found to exist between the coefficients r 2 and r 3 , namely: r 3 = −ir 2 k 0z /k x . From the magnetic field components above, we can use Eq. (7b) to easily deduce the electric field components for the vacuum region: where For the WS region, when implementing Maxwell's equations, the anisotropic nature of the WS must be explicitly taken into account by involving Eq. (1) in the calculations. The general solution to the E field in the WS region should then be a linear combination of four wavevector components To determine the set of coefficients {r 1 , r 2 , r 3 }, and {a 1 , a 2 , a 3 , a 4 }, it is necessary to invoke matching interface conditions and boundary conditions. But first, the remaining E and H fields must be constructed. This can be accomplished by first using Eq. (7b) to write the following relations: While the other Maxwell equation, Eq. (7a), gives, where the translational invariance in the x-direction has been used for the derivatives: Upon matching the tangential electric and magnetic fields at the vacuum/WS interface, and using the boundary conditions of vanishing tangential electric fields at the ground plane, it is straightforward to determine the characteristic equation that governs the EM modes of our structure. The resultant transcendental equation that must be satisfied can be written as D = 0, where, Here we have introduced the dimensionless quantities: For the lossless guided wave modes studied in this paper, the dispersion equation (Eq. (18)) is either purely imaginary or real. In the absence of gyrotropy (γ = 0), Eq. (18) becomes decoupled, and the corresponding dispersion equation for purely diagonally anisotropic media reduces to, 53 and, The allowed modes in Eq. (19) (20) accurately account for the guided mode solutions. When the gyrotropic parameter is negligible, an incident EM wave interacting with the WS retains its initial polarization state, thus limiting complications from Faraday and Kerr effects. Of course, as discussed above, when γ is not negligible, all three polarization directions must be taken into account, and nontrivial coupling can arise between the components of the E and H fields.
An insightful quantity that yields the fraction of EM energy that is confined to the WS region is the confinement factor η, defined as, where U i ≡ Ωi dz U i (z) is the energy density U i integrated over region Ω i (i = 0, 1). The energy density, which accounts for any possible anisotropy and dissipation present, is written, which ensures that the energy density is positive, as required by causality. Related to the energy density, is the time-averaged energy flow, given by the real part of the Poynting vector S, which for energy flow along the direction of the interface is expressed as, The corresponding power in each region P xi is then defined in the usual way as an integral of the Poynting vector over region i: P xi = Ωi dzS xi . From the Poynting vector and energy density defined above, it is possible to construct another physically relevant quantity, namely the energy transport velocity, 60,61 v T , which is the velocity at which EM energy is transported through a given region of the waveguide structure, and defined as: where the brackets indicate spatial averaging over all space, and only time-averaged quantities are considered here.

C. Discussions and numerical results
Having established the methods for determining the waveguide modes for the WS structure, we now consider a range of material and geometrical parameters that leads to the propagation of localized EM waves along the WS surface. We consider subwavelength widths d/λ 1, leading to a decoupling of the EM field components and excitation of TM modes where the components E x , E z , and H y dominate. We search for modes that are in the parameter regime where the dissipative components of zz and vanish. From the previous discussion in Sec. II A, this implies that we focus on frequency bands that satisfy ω < 2µ/(1+β). The process for calculating the modes is as follows: For a given frequency ω, and relevant WS parameters, each component of the tensor is calculated. These components are inserted into the transcendental equation (18) to determine the permitted propagation constants κ x that are associated with the chosen ω. This creates a map of a finite number of dispersion bands for a given waveguide width d that ensure the boundary conditions are fully satisfied. Finally, from the dispersion diagram, the EM fields, and energy characteristics can be straightforwardly extracted. In presenting the results, we generally scale ω by the energy unit v F |Q|, and define ω ≡ ω/(v F |Q|). Similarly, when considering the chemical potential in the WS, we take the representative dimensionless value corresponding to µ ≡ µ/(v F |Q|) = 0.2. We also set the energy cutoff to Γ ≡ Γ/(v F |Q|) ∼ 8, so that the linearized model remains applicable.
In Fig. 2(a), we present the dispersion diagram in terms of the dimensionless propagation constant κ x and frequency ω. Various points are numerically labeled on the diagram for discussion below. We consider subwavelength widths corresponding to k c d = 0.03, where k c is a characteristic wavenumber corresponding to k c = ω c /c, at the normalized design frequency ω c = 0.12. The conical tilt is fixed at β = 0.9. As seen, there are a discrete number of bands, which arise from the geometrical constraint placed upon the fields from the bounding surfaces. Also, since the fields in the vacuum decay over the length scale k −1 0 , regions of the dispersion curves near the light line (κ x ∼ 1) correspond to a larger decay length, and thus the fields reside mainly outside the guiding layer. The curves that lie to the left of the vertical dashed line identify modal solutions that are predominately evanescent, with the EM field decaying on both sides of the WS/vacuum interface. This arises when both and zz are negative. For frequencies to the right of the vertical line, the modes transition from predominately evanescent in the WS, to sinusoidal, whereby > 0 and zz < 0. The modes at these frequencies of course still decay in the vacuum re-gion. In Fig. 2(b), the normalized wavevector in the WS, κ 1z , is shown. Due to the decoupling of the TE and TM modes, one wavevector is dominant, and so we have κ ± → κ 1z . At the transition point between oscillatory and evanescent modes (vertical line at ω ≈ 0.117), ≈ 0, and consequently the wavevector drastically diminishes, corresponding to a very long wavelength in the WS. On the other hand, when zz ≈ 0 (point 3), the wavevector is noticeably larger. The correlations between the mode diagrams and the EM response is shown in Fig. 2(c). Notably, and zz become vanishingly small at different frequencies, leading to an interesting anisotropic ENZ response. In particular, at the crossover point ω ≈ 0.117, ≈ 0, and zz < 0. Indeed, when the energies associated with the chemical potential and frequency are similar, ≈ 0 when ω ≈ 2µ/ 3π/α + log(4Γ 2 /(3µ 2 )). The other ENZ scenario occurs at the higher frequency ω ≈ 0.18, where zz ≈ 0, and > 0. This is also in the vicinity where the slope of the dispersion curve in Fig. 2(a) changes sign (labeled 3). This feature of the dispersion diagram corresponds to where κ x close to the light line bends back and rapidly increases. Similar behavior is seen in general, when a waveguide is in contact with another medium possessing permittivity components that are opposite in sign. The corresponding energy flow along the interface that separates the two media undergoes an abrupt reversal when going from vacuum to the WS or vice versa. Additionally, extreme field enhancement can occur when either of the permittivity components is near the ENZ regime. For example, since the boundary conditions dictate that the ratio of the normal components of the electric fields at the semimetal interface are E z1 /E z0 ∼ 1/ zz , a substantial mismatch in field strengths can clearly occur when zz ∼ 0. Although its effects are weak for a subwavelength WS, for completeness the inset presents the frequency dependence to the off-diagonal component γ calculated via the methods described in the Appendix. Lastly, in Fig. 2(d), the energy confinement factor η is shown as a function of normalized frequency. The labeled numbers identify the corresponding (κ x , ω) pairs in Fig. 2(a). The fraction of timeaveraged EM energy within the WS at frequency points 1 and 2 is shown to be extremely small due to the fact that modes with κ x ∼ 1 correspond to the EM fields residing mainly in the vacuum region. The upper branch at points 3 → 5 reveals that a substantial percentage of the energy density is contained within the semimetal. We see that η is maximal at position 3, where the slope of the dispersion diagram changes sign [see Fig. 2(a)], and zz ≈ 0, giving rise to a large increase of the electric field normal to the interface.
To gain additional perspectives on the energy flow in the system, we present in Fig. 3(a), the energy transport velocity v T as a function of the normalized frequency. To calculate this quantity we take the ratio of the Poynting vector to the energy density per Eq. (24). For the subwavelength structure under consideration, the time-averaged energy density can be calculated approximately for the WS region by: For the vacuum region, we have: Likewise, the time-averaged Poynting vector component S x for d/λ 1 can be expressed as, and for the vacuum region: In the absence of loss, the energy transport direction lies solely along the interface (x-direction), so that v T = v Tx . Note that when κ x → 1, v T corresponds to the expected phase velocity, or velocity at which plane wavefronts travel along the xdirection: v T = c/κ x . It is evident from Fig. 3(a), that v T obeys casualty, with v T ≤ c. Due in part to the competing effects that can arise in structures with counter-propagating energy flows on opposite sides of the waveguide interface, the The total normalized power of the system. Both vT and Ptot are calculated for the modes described by the first dispersion curve in Fig. 2(a), and correlated with the number labels. The WS has its conical tilting parameter set to β = 0.9.
transport of energy can occur at effective speeds smaller than c. In particular, it is evident that v T can become very small for the modes labeled 3 → 5. Note that from the inset, we see that over the entire frequency range, v T is slightly negative, but never zero until at the crossover point (3), where the dispersion curve in Fig. 2(a) changes direction. For modes 1 and 2, where κ x ≈ 1, a large portion of the EM field resides in the vacuum region, and thus the net transport of energy has v T /c ≈ 1. Slow-light phenomena have been explored previously 62,63 in ENZ and metamaterial systems. The flow of energy involves both the interplay between the propagating and stored forms of energy. The proposed WS waveguide thus serves as an effective platform for studying slow-energy effects in dissipationless systems, where the energy propagation in adjacent regions are oppositely directed due to the sign change of the permittivity across the interface. Having tunable control over the speed of energy transport can also have practical uses in optical memory devices, biosensing 65 , and chemical sensors. 64 As a complimentary study, we show in Fig. 3(b), the total power P tot of the system, where we define, The upper branch of the power curve at unity (modes 1 → 2) corresponds to the dispersion curve where κ x ∼ 1 in Fig. 2(a), where the energy flows predominately outside the semimetal. The net power drops drastically to zero at the transition point (mode 3), where the slope of the dispersion curve changes sign. This follows from the superposition of two counterpropagating energy flows, and is consistent with the behavior of the group velocity v g ≡ ∂ω/∂k x , which becomes vanishingly small at point 3. Reducing the frequency, and going from points 3 → 5, it is apparent that the power reverses and now has a net flow in the opposite (negative) direction, corresponding to the upper branch of the dispersion curve in Fig. 2(a), until the net power again vanishes at ω ≈ 0.112. As alluded to earlier, for frequencies in which the diagonal permittivity components change sign when crossing into another medium, the resultant flow of energy in each region can often times be oppositely directed. To visualize the directional dependence to the energy flow, we present in Fig. 4, the spatial behavior of the Poynting vector in the surrounding vacuum region (z/d > 0), and within the WS (0 ≤ z/d ≤ −1). We consider the frequency and propagation constant corresponding to point 3 on the dispersion curve in Fig. 2(a). At this point zz is slightly negative, with zz ≈ 0, and ≈ 1. First, in Fig. 4(a), the instantaneous Poynting vector is shown. We see that the energy flow tends to form a vortex-like pattern in which the flow of energy upward is counteracted by the flow downward. Since the normal component of the E field is discontinuous at the interface, the energy flow undergoes an abrupt change at z = 0. Since there is no dissipation however, when averaging over a complete period, there is no net power flow in the z direction. This is seen in Fig. 2(b), where the time-averaged Poynting vector is presented. The cancellation of the energy flow normal to the interface leads to a net propagation of energy only along the interface (x direction). We saw earlier in Fig. 3(b) that this mode corresponds to a vanishing net power flow. Thus, the spatially integrated Poynting vector component along the interface in each region completely cancel one  Fig. 3(a) where the energy velocity vanishes (labeled 3). The interface between the vacuum region and the WS is located at z = 0.
another, leading to a net effective energy velocity of zero [see Fig. 3(a)].
We now turn our attention to the explicit behavior of the fields themselves. In Fig. 5(a) we show the spatial behavior of the electric field parallel to the interface E x . Four representative frequencies are considered along the dispersion curve show in the inset. These modes correspond to the regime that is left of the vertical dashed line in Fig. 2(a), and thus these modes are strictly evanescent on both sides of the interface. It is noted that as the frequency approaches the lowest frequency case, ω ≈ 0.112 (labeled 1), the electric field becomes strongly enhanced. As was observed in Fig. 3, this correlates with the frequency in which v T ≈ 0, and the slope of the dispersion curve vanishes. Note that over the considered frequency range, both diagonal components of are negative, and ≈ 0. Thus, these are purely evanescent modes with field profiles that decay on both sides of the vacuum/WS interface, as is clearly seen in Fig. 2(a). The energy density, which contains both the electric and magnetic fields also exhibits strong field localization at the interface as seen Fig. 2(b). Next in Figs. 5(c) and 5(d), color maps illustrate the full spatial profiles the tangential (E x ) and normal (E z ) electric field components, respectively. Here the propagating term e ikxx is included for all field components. The tilt parameter β is set to β = 0.9, while ω is fixed to the frequency labeled 1 in the mode diagram shown in the inset of Fig. 5(a). It is evident from Figs. 5(c) and 5(d) that the propagating electric field has both field polarizations strongly confined to the interface. As shown in Figs. 5(a) and 5(b), the EM fields also become sub- stantially enhanced at this frequency.
We now discuss how the tilt of the Weyl cones directly influences the mode diagrams and corresponding EM field characteristics. In Fig. 6(a), the dimensionless propagation constant κ x is shown as a function of tilt β. For clarity of discussion later, the numbers label relevant parts of the dispersion diagram. In Fig. 6(b) the wavevector in the WS, κ 1z is also shown over a wide range of β. In contrast to what was observed in Fig. 2(b) when varying the frequency, here the wavevector in the WS monotonically increases along the dispersion path 1 → 5. In Fig. 6(c) the EM response of both the diagonal (main plot) and off-diagonal (inset) components of is shown over the same range of β. For the given frequency, ω = 0.117, both components and zz are negative, with ≈ 0. Examining Fig. 6(a), we see that the main dispersion curve has the lower branch with points 1 → 3 close to the light line. The upper branch from 3 → 5 is associated with substantially larger propagation constants that increase rapidly with tilt. The corresponding confinement factor is presented in Fig. 6(c). As expected, the fraction of energy η contained in the WS is very small for points 1 and 2, where κ x → 1, reflecting the predominance of the energy density outside of the guiding semimetal layer. There is a significant fraction of the energy density in the WS for modes occupying the upper dispersion curve along the path 3 → 5. Moreover, η → 1 for conical tilts in the vicinity of points 3 and 5, where ∂β/∂κ x → 0 in Fig. 6(a). Thus, we see that the The dominant wavevector in the WS, κ1z, is shown in (b). The electromagnetic response of the WS waveguide is shown in (c), where the diagonal components of the permittivity tensor and zz are presented in the main plot, while the inset reveals the off-diagonal component γ. The fraction of EM energy η that is confined to the semimetal for the modes calculated in (a) is shown in (d). The normalized frequency corresponds to ω = 0.117, which is the point labeled 5 of the dispersion curve shown in Fig. 2(a). Note that κx and κ1z are shown on a logarithmic scale.  Fig. 6(a) (with the labels 3 → 5). The frequency is identical to that used in Fig. 6. tilt of the Weyl cones plays an important role in determining the propagation and spatial localization properties for guided wave modes.
To study the guided modes in the ENZ regime further, we keep the wave frequency the same as in Fig. 6, corresponding to where ≈ 0, and investigate the spatial behavior of the fields. Thus, in Fig. 7, the components of the electric field E x and E z are shown as a function of the scaled coordinate z/d for varying tilt of the Weyl cones. The three values of the tilt parameter shown in the legend corresponds to the points 3 → 5 of the mode diagram in Fig. 6(a). Notably, there is a linear variation in E x , while E z is constant throughout the WS layer. The E-field for both components then decays once en-tering the vacuum region, as is required for guided waves. For both components, increasing the tilt results in a field enhancement, as zz becomes increasingly more negative, resulting in modes with larger κ x [see upper branch of Fig. 6(a)]. To understand this overall behavior explicitly, we write the general expressions for the EM fields and then take the limit → 0, arriving at: where the factor e ikxx has been suppressed. The linear dependence seen in Fig. 7 for E x1 is clearly evident from Eq. (30), while the spatial dependence in E z1 and H y1 drops out, leading to the uniform component of the electric field normal to the interface, seen in Fig. 7(b). Interestingly, despite the fact that E z1 is a function of κ x and zz , both of which depend strongly on β [see Fig. 6(a) and 6(c)], they tend to counteract one another, leading to the weak β dependence observed for E z1 in Fig. 7(b). It is instructive to also examine the energy flow via the Poynting vector. Within the ENZ regime, we find, which is independent of position in the WS. The spatially averaged energy density in the WS is, To obtain the speed of EM energy transport, it is then a simple matter to take the corresponding ratio via Eq. (24). The dispersion diagram can also be found analytically in the ENZ regime. When ≈ 0, the propagation constant can be writ-ten: which is a convenient expression that negates the need to numerically solve a transcendental equation for the dispersion diagram.

III. CONCLUSIONS
We have analyzed the dispersion diagrams and electromagnetic field behavior of a subwavelength planar Weyl semimetal waveguide system where the Dirac cones in the touching points of valance and conduction bands possess the possibility to be tilted. We found that the electromagnetic response for each diagonal component of the permittivity tensor can be made is vanishingly small, creating an anisotropic epsilon-near-zero configuration. The admitted waveguide modes were presented in the form of dispersion diagrams, where regimes of high energy-density confinement in the Weyl semimetal layer were found. We also showed how the net system power can be manipulated to vanish at critical points of the dispersion curves, where the propagation velocity of electromagnetic energy is zero. Our results offer waveguides made of Weyl semimetals with controllable energy flow velocity and direction by materials and geometrical parameters.
and γ (s) where the following limit is taken in the integrals: ω k → ω + iδ. This ensures proper convergence when numerically evaluating Eqs. (A1)-(A2). We have also defined ζ s,t ≡ tp + p z β s , p = p 2 z + p 2 ⊥ , and a momentum cutoff along the z axis, Γ. Generally, the cut-off Γ is a function of the tilt parameter. Nevertheless, in our calculations, we choose a large enough cut-off and neglect the contribution of β to Γ. The cut-off Γ 0 > v F |Q| is introduced for the correct definition of the vacuum contribution.