On the possibility of superluminal energy propagation in a hyperbolic metamaterial of metal-dielectric layers

Metamaterials usually refer to artificially designed structures consisting of periodically arranged resonators,¹ dielectric-metal layers² or intrinsic resonant components for achieving unusual electromagnetic properties such as negative refractive index or indefinite permittivity.^3,4 With such unusual properties, novel devices such as superlens⁵ and hyperlens⁶ for subwavelength imaging and invisibility cloak⁷ for hiding objects can be realized experimentally in practice. The resonance and directional-conductive characters of these metamaterials imply that they are inherently dispersive, absorptive, and usually anisotropic. A great many research works have already been done for understanding the propagating properties of electromagnetic waves in these media, and various novel devices have been proposed and tested.^8–12 However, studies concerning how wave packets or signals propagating in these dispersive media are few.¹³ The wave packet propagation problem is theoretically interesting because it is related to some fundamental problems like the causality or superluminality. It may also be technically relevant because with the knowledge about how wave packets propagate and distort in the media, we may develop useful techniques for controlling or manipulating the signals evolving in these media, which is highly desirable for applications related to time dependent phenomena.

It is well known that any kind of wave propagation process must satisfy causality.¹⁴ That is, the cause must happen before the appearance of the result. Special relativity tells us that if superluminal signals can be sent as messages connecting event A and event B in an inertial frame, then it is always possible to find an inertial frame in which the observed orders of time for them are exchanged, violating the constraint of causality. Another restriction deduced from relativity is that no massive particle can move faster than the speed of light in a vacuum.¹⁵ These considerations imply that waves can never have superluminal propagation velocity for the energy or information carried by them. However, recent studies concerning the propagation of electromagnetic,¹⁶ acoustic,¹⁷ and quantum waves¹⁸ in various kinds of dispersive and dissipative media revealed that, under suitable conditions and through various mechanisms, superluminal phenomena do exist and can be observed, though in some subtle ways they never break causality and do not violate the restriction of speed limitation from relativity.¹⁹ 

Up to now, most studies about the superluminal phenomena considered isotropic media or one-dimensional propagation.^20–26 About a decade ago, a study claimed that the group velocity of hyperbolic propagating mode can become superluminal under suitable conditions.²⁷ However, a recent theoretical work concerning passive linear media with negligible loss proved that the group/energy velocity of the electromagnetic propagating modes can never exceed the speed of light in a vacuum.²⁸ These contradictive results attracted our attention and inspired us to further explore the propagation behaviors of wave energy in a dispersive-anisotropic medium.

In this paper, we provide a rigorous proof showing that the group velocity of any propagating mode inside the effective medium that derived from the periodic structure of dielectric and metallic layers, no matter whether it obeys elliptical or hyperbolic dispersion, is always slower than the speed of light in a vacuum. This is the most important result of this paper. We also numerically simulate the wave packet propagating behavior based on the Fourier transform method and the energy density formula derived in this paper. Although our result on the mode propagating velocity confirmed the conclusion given in Ref. 28, our proof is concrete in details and directly based on the material properties, different from that of the more abstract approach in Ref. 28. Besides, though in the lossless limit our time dependent energy density formula applying to a propagating mode is the same as that derived from the corresponding formula in Ref. 28, for the case containing nonzero loss our energy density formula follows our previous work,²⁹ excluding the “dissipation part” in the formula of Ref. 28 that represents the time integration of the power loss. In addition, very detailed results about the signatures and dispersions of the principal values of permittivity tensor are given, which is rare in existing publications.

The metamaterials discussed in this paper can be realized effectively in the long-wavelength limit by the one-dimensional photonic crystals consisting of alternatively arranged dielectric and metallic layers. In Sec. II A, we start with a one-dimensional photonic crystal (1DPC) of dielectric-metal layers and show that the principal values of the effective permittivity tensor along the two directions parallel to the material layers is of the same form as in the Drude model, implying a plasma behavior, whereas in the direction perpendicular to the layers, the effective medium mimics a Lorentz medium. We then derive the formula of group velocity for propagating mode in Sec. II B and discuss in Sec. II C the restrictions on its propagation directions and prove rigorously in Sec. II D that they must be slower than the light speed in vacuum if the dielectric materials used in this 1DPC is a real material having relative permittivity larger than one. We also derive the formula of the dynamical wave energy density in Sec. II E for arbitrary time varying electromagnetic fields in the effective medium. The energy formula is then used in the numerical study of the energy propagation velocity of a wave packet in Sec. III A and III B. The wave packet is formed by a supposition of many plane waves with properly chosen coefficients for different propagation directions and frequencies. We then compare the numerical results for the wave packet and the theoretical results for the propagating mode and discuss the reasons of their differences. The conclusion is presented in Sec. IV.

To begin with, we consider a one-dimensional periodic structure which has one metal layer and one dielectric layer per unit cell, as is shown in Fig. 1. In the long wavelength limit this photonic crystal behaves like an effective anisotropic medium which in some frequency ranges has indefinite eigenvalues of its permittivity tensor. We arrange these layers to be parallel to the xy plane and periodic along the z axis, so if we define xz plane to be the plane of incidence then the electromagnetic (EM) fields in this structure can be set to be independent of the y coordinate. Since the anisotropic feature appears only for the TM (E_x E_z H_y) polarized waves, hereafter we consider only waves of this polarization.

Denoting the filling fraction of the dielectric layer in one unit cell as F, and ε_d and ε_m the permittivity of the dielectric and metallic layer respectively, the principal values of the effective permittivity tensor of the metamaterial along the x and z directions under the long wavelength limit are given by:

We assume that the permittivity ε_d of the dielectric layers is a constant and the metal permittivity ε_m is described by the Drude model with negligible dissipation effect, that is:

where ε₀ is the dielectric constant of the empty space.

The symbol ω_p stands for the plasma frequency which depends on the mass and concentration (number density) of the electrons in the metal. Substituting Eq. (2) to Eq. (1), we get

where $ε̃x$ and $ε̃z$ are the principal values of the relative permittivity along the x and z axes. The coefficients A_x and A_z are given in Eq. (4). We observe that under the long wavelength limit, the effective permittivity in the x direction still has the Dude-model form, whereas along the z direction the metamaterial behaves like a polaritonic medium having Lorentz-type dispersion. That is, along the z direction there is a frequency gap between two resonance frequencies ω₀ and ω_p. Hereafter we denote the relative permittivity of the dielectric layers as $ε̃d=εd/ε0$⁠. The parameters in Eq. (3) are given by:

Here both $ωpeff$ and ω₀ are functions of F, and they become the same if F = 1/2.

In this section we derive the formula for the group velocity of a propagating mode in a general anisotropic dispersive medium. This velocity can be further shown to be equal to the energy velocity defined by the ratio between the time averaged Poynting vector and the time averaged energy density. We will not derive this result here but only recognize it as a fact.

The dispersion relation for TM propagating mode in an anisotropic dispersive medium characterized by the principal values of permittivity $εz=ε0ε̃z$ and $εx=ε0ε̃x$ can be derived from Maxwell’s equations as

The group velocity v_g = ∇_kω for this mode can be deduced by taking the k-gradient of the dispersion relation:

Since the principal values of the permittivity of the medium are frequency dependent, so there are the relations $∇kε̃x=∂ε̃x∂ω∇kω$ and $∇kε̃z=∂ε̃z∂ω∇kω$⁠. Substituting them into Eq. (6), we get

Multiplying Eq. (7) by ω and using the dispersion relation (5), we have

thus the group velocity in the dispersive anisotropic medium is:

In this section we discuss the existence condition and the restrictions on the propagation direction of a propagating mode. To be more specific, we define

For a propagating mode all of the three quantities must be real, and conventionally k takes positive value. Note that according to Eq. (5) if both $ε̃x$ and $ε̃z$ were positive there would be no restriction on the propagating direction. On the other hand, if both of them were negative there would be no propagating mode at all. For the interested indefinite situation $ε̃xε̃z<0$⁠, we have either $ε̃x<0$ or $ε̃z<0$⁠. Substituting Eq. (10) into Eq. (5), we get

Suppose we have $ε̃x>0>ε̃z$⁠, then Eq. (11) implies that the propagating mode can only have wave vector oriented inside the angular region $|θ|<sin−111+|ε̃x/ε̃z|$⁠. Similarly, if we have $ε̃z>0>ε̃x$⁠, then the wave vector must oriented in the region $|θ|>sin−111+|ε̃x/ε̃z|$⁠.

The signs of the principal values of permittivity defined in Eq. (3) are controlled by the wave frequency ω and the values of ω₀, ω_p, and $ωpeff$⁠, and these characteristic frequencies are determined by the filling fraction F through Eq. (4). For simplicity we define “boundary angle” as $θb=sin−111+|ε̃x/ε̃z|$ and assume $ε̃d=εd/ε0>1$⁠. For F belonging to different ranges of frequency we have the following three tables:

These tables summarize the existence conditions for propagating modes and the restrictions on their propagating directions. The F = 0 and F = 1 cases respectively correspond to pure metallic and pure dielectric materials. In the next section we will show that if $ε̃d=εd/ε0>1$ were assumed, any propagating mode can only have group velocity slower than the speed of light in vacuum.

To discuss the speed of a propagating mode, we define a factor β_g = |v_g|/c, which represents the ratio between the group speed and the universal speed of light in vacuum. If β_g can be proved to be always smaller than 1, then the subluminality is established. According to Eq. (9) and Eq. (10), β_g takes the value:

Which can be rewritten as $βg=γη$⁠, where

If both γ and η can be proved to be smaller than 1, than subluminality must be true. Using Eq. (3), we find:

Note that according to Eq. (4), if we assume $ε̃d=εd/ε0>1$ then we have A_x > 1 and A_z > 1, and therefore the two factors in Eq. (14) are also larger than 1. This result directly leads to γ < 1 because its denominator is larger than the numerator.

For the η factor, there are three situations we need to consider. First, if both $ε̃x$ and $ε̃z$ are positive, we have $ω2ε̃x′2ω>ε̃x$ and $ω2ε̃z′2ω>ε̃z$⁠, and thus η < 1. In the second situation, we have $ε̃x>0>ε̃z$ and a positive numerator of η. Now if we replace the original numerator with the positive term $ε̃xkzε̃x2$⁠, we still get a ratio smaller than one. Therefore we know in this situation we also have η < 1. The third situation is $ε̃z>0>ε̃x$⁠. We can just exchange the role of $ε̃x$ and $ε̃z$ in the second situation and follow the same reasoning to get η < 1. There is no need to consider the case that both $ε̃x$ and $ε̃z$ are negative, because under this situation there is no propagating mode. Since for all the situations discussed above we always have γ < 1 and η < 1, we conclude that $βg=γη<1$ and the subluminality for any propagating mode is established. Note that in the above discussion we always assume the condition $ε̃d=εd/ε0>1$⁠. Without this assumption, we would not be able to establish the subluminality in general. We thus get the main result of this paper: for the effective medium derived from the one-dimensional photonic crystal consisting of dielectric-metallic unit cell, if the dielectric material is a real material having relative permittivity larger than 1, then the mode velocity in this effective medium must be slower than the speed of light in vacuum. Numerical results for filling fraction F = 0.25, 0.4, 0.75 and 0.95 are shown in Fig. 2. Here the “k-angle” denotes the angle (in degrees) between the wave vector k of the propagating mode and the z-axis (the direction of the periodicity), and the color in color bar represents the group velocity to light speed ratio β_g defined in Eq. (12). The blank regions represent the forbidden combinations of (ω,k) for propagating modes. As one can see, all the results confirm the conclusion of subluminality.

Based on this discussion, now it’s easy to understand that the superluminality of the mode or wave packet velocity in Ref. 27 is caused by the fact that the frequency dependence of the permittivity tensor were not taken into account in the derivation of the group velocity formula, equation (7).

In the previous subsection, we have rigorously proved that the group velocity of a propagating mode in the effective hyperbolic metamaterial we defined in II A must be subluminal. It is well known that for a propagating mode of frequency ω one can define its energy velocity as $ve=SU$ and this velocity can be proved to be equal to the group velocity, i.e, v_e = v_g if the absorption effect can be neglected. Here $⟨A⟩$ stands for the time average of the time-dependent quantity A(t), S = E × H is the Poynting vector, and U is the energy density of the electromagnetic wave in the medium. However, instead of proving this relation, in the following sections we will compare v_g with the traveling speed of a wave packet which has central frequency ω and formed by a weighted sum of many propagating modes. We thus derive the instantaneous energy density formula here and apply it in the following sections for studying the energy propagation problem.

In classical electrodynamics the displacement field D(r,t) is related to the electric field E(r,t) and electric polarization field P(r,t) through the relation:

For the purpose of deriving the energy density formula of the EM field in the present medium, we find this relation is not very convenient to use. Instead of using Eq. (15) directly, we define a Q(r,t) field through the relations:

Now we derive the dynamical relations between the Q(r,t) field and the applied electric field. Multiply both sides of Eq. (3) with −ω²E_x and −ω²E_z, apply Eq. (16) and replace −ω² with ∂²/∂t², we get:

In fact, the terms A_xε₀E_x and A_zε₀E_z in Eq. (16) represent the D_x and D_z at infinite frequency, as can be checked from Eq. (3). They correspond to the non-dispersive response of the effective medium to the applied electric field. On the other hand, the Q(r,t) field represents the dynamical polarization field (the dispersive P field) of the medium when applying a time-varying electric field. Note that in this section all the quantities mentioned such as the electric field or the Q-fields etc. are real value fields, not phasors (complex quantities).

Substituting Eq. (17) into Maxwell’s equations, a Poynting theorem equation can be obtained. The left-hand side of this equation is the negative divergence of the Poynting vector −∇⋅S = −∇⋅(E × H), whereas the right-hand side is $E⋅∂D∂t+H⋅∂B∂t$⁠. Similar to the previous works,^29,30 the right hand side of the equation can be identified as the time rate of change of the field energy plus the power loss. Since in this paper we assume the power loss of the medium is small enough to be negligible, the right-hand side can be written as $∂U∂t=∂Ue∂t+∂Ub∂t$⁠. Here $∂Ue∂t=E⋅∂D∂t$ is the rate of change of the electric energy density including both the contributions from the electric field itself and the material response of the medium. Similarly, $∂Ub∂t=H⋅∂B∂t=∂∂t12μ0H2$ is the rate of change of the magnetic field energy density. Finally, U represents the total energy density to be identified. It is obvious that the magnetic energy density is simply given by $Ub=12μ0H2$⁠, but the electric part U_e is not that simple and need to be identified.

According to Eq. (17), we have:

thus the electric field energy density with material response should be identified as:

The electric energy density U_e in Eq. (19) can be divided into two parts which can be recognized as the contributions from the electric field and the non-dispersive response of the medium (the first two terms on the right hand side of Eq. (19)) and the dispersive responses of the matter (the other terms). Once we know both the E and H fields in the medium, we can calculate U_e from Eq. (19) and $Ub=12μ0H2$⁠, thus the total energy density U = U_e + U_b is obtained. The dispersive contribution (the second and third terms on the right hand side of Eq. (19)) of the energy density are from the kinetic and potential energies of the vibrating dipoles in the medium. The dynamical vibration of the dipoles is induced by the electric field in the medium.

In the next section we will numerically study the propagation of a wave packet in the hyperbolic medium. Two numerical velocities will be compared: the wave-amplitude velocity and the wave-energy velocity. The former is defined by tracking the motion of the amplitude center of the wave packet, whereas the latter is defined by tracking the motion of the energy-center of the wave packet, based on the formula derived in this subsection.

We can build a wave packet of arbitrary shape by adding the plane waves of propagating modes with appropriate weight function of frequency and wave vector. This is relying on the fact that the plane waves e^{i(k⋅r−ωt)}with k and ω satisfying the dispersion relation Eq.(5) form a complete set, so any solution of Maxwell’s equations in the medium can be written as a superposition of them. For numerical simulation of TM wave packet, we sum the magnetic fields and assume that the weight function has the following Gaussian form:

The ω_b and θ_b determine the distribution width of the plane wave bases in frequency and in propagating direction, respectively, whereas the ω_c and θ_c represent the central frequency and the central propagating direction (hereafter we call it incidence angle) of the wave vectors. The wave packet consisting of m frequencies and n wave vectors is thus written as:

In practice, the number of wave bases for constructing the wave packet should be large enough to prevent the possibility of more than one packets appear in our interest spatial region of observation.

As mentioned before, the signs of $ε̃x$ and $ε̃z$ are determined by the wave frequency ω and three characteristic frequencies ω₀, ω_p, and $ωpeff$⁠, and these parameters are determined by the filling fraction F through Eq. (4), as summarized in tables A to C. According to these results, the propagating directions of the modes are restricted by the condition |θ| > θ_b or |θ| < θ_b, or without any restriction. To simplify our discussion without reducing relevant generality, we assume $F=12$⁠, i.e., the case in Table C. Hereafter we also assume $ε̃d=2$⁠. Under this choice the $ε̃x$ and $ε̃z$ as functions of frequency are shown in Fig. 3.

Notice that in Fig. 3, there are three different frequency bands characterized by the signs of the permittivity pair $(ε̃x,ε̃z)$⁠: (−, +), (+, −), and (+, +). Similarly, the theoretical group velocity (Eq.(12)) as function of the phase propagating angle (the wave vector direction of the mode) and frequency also exhibits three different patterns, as shown in Fig. 4. The similar behaviors of the permittivity tensor and group velocity reflect the dispersive and anisotropic features of the effective medium.

In this research, we are not interested in the third region of Fig. 3 because in it both of the two principal values the permittivity are positive. We will only consider the case of indefinite medium with hyperbolic dispersion relation. Since the first and second regions are not very different, we can choose either one of them. Hereafter we choose the second region as the frequency domain and do all the numerical simulations. Under this choice the effective medium has positive ε_x and negative ε_z.

Now we can build the wave packet using Eq. (21). We define the wave packet center just like one defines the center of mass in a dynamical problem, but for the present problem the mass distribution is replaced by the magnitude of the instantaneous H-field or the energy density of the electromagnetic fields (including the response energy of the medium). We trace the motion of the center of the wave packet numerically, and the numerical “group velocity” of the wave packet is defined as the group center displacement in the time interval of one step of time counting. The typical field energy density distribution at three different times is shown in Fig. 5. In this simulation we choose ω_c = 1.15 GHz, ω_b = 0.14 GHz, θ_b = 20°, m = n = 25(625 plane waves), and the energy density distribution for t = 0, 15, 30μs are shown in the figure.

The comparison between the mode group velocity (Eq. (9)) and the numerical wave packet velocities in the indefinite frequency domain for central phase propagating angle of 10 degrees and 30 degrees are shown as Fig 6. From it we find that for the wave-packet propagation, no matter whether we evaluate its velocity based on the amplitude or energy distribution, the result is very close to the theoretical group velocity of the mode if the bandwidth of the wave packet is not too wide. However, if the bandwidth of the wave packet is a little wider, the two possible packet velocities deviate a little bit from the theoretical group velocity of the mode. This is because some Fourier components the wave packet become evanescent waves and cannot penetrate through the media, and this leads to the shift of the central frequency of the wave packet.

In Fig. 7, we compare the velocity diagrams of the theoretical group velocity of the propagating modes and the numerical wave packet velocity defined with the energy density scheme. In constructing the numerical velocity diagram, various cases including 30 central frequencies and 20 central phase propagating angles are considered. The similarity between the theoretical and numerical results is very obvious.

In this paper, the effective electromagnetic properties in the long wavelength limit of the 1DPC of dielectric-metal layers is studied. This metamaterial behaves like a plasma or Lorentz medium with respect to the applied electric fields parallel or perpendicular to the layers, respectively. Based on the derivation of the Poynting theorem, the formula for the dynamical energy density is obtained. This dynamical energy includes two parts: the energies from the electric and magnetic fields themselves, and the kinetic plus the potential energies from the charge carriers in the medium. For the propagating modes in the medium, we derived the restrictions on the propagating directions of them. According to our theoretical analysis and numerical results, the group velocity of a propagating mode in the effective medium is always subluminal (slower than the speed of light) if the Drude and Lorentz dispersions already lead to subluminal propagation separately, that is, they are real materials. Thus our results confirm the conclusion of the theoretical work in Ref. 28 and reveal the reason of why superluminal propagation were derived in Ref. 27. Besides, for the wave-packet velocity, no matter how we evaluate it, based on tracking the packet amplitude or energy distribution, the result is very close to the theoretical value of the group velocity for the propagating mode having the frequency close to the central frequency of the packet if the bandwidth of the packet is not too wide. If the bandwidth of the wave packet is a little wider, however, the two possible wave packet velocities (amplitude and energy) become deviated from the theoretical group velocity of the mode. We believe this is because some Fourier components of the wave packet become evanescent and cannot penetrate through the media, and this leads to the shift of the central frequency of the wave packet in the propagation process.

Our conclusion that no superluminal propagation of wave energy can happen in a dispersive-hyperbolic metamaterial is based on two main assumptions. The first is that the dielectric layers must have permittivity greater than 1, and the second is that the metal layers have permittivity of Drude-model type dispersion. Both of which are typical dielectric properties for real dielectric and metallic materials, respectively. It is in this sense we call them “real materials” in our paper. In addition, we assume through this paper that the medium is lossless. This assumption may seem unphysical because the Kramers-Kronig relations require a nonzero positive damping coefficient in the permittivity formula.³¹ However, this requirement does not give a lower bound to such damping coefficient but only restricts the sign of it. Thus a lossless medium is not unphysical because it can be treated as an absorptive medium with a negligible damping coefficient. Finally, we want to stress that our conclusion does not exclude some superluminal phenomena happen in finite thickness gain medium or evanescent regions.³² Our results are about wave propagation in an infinite or semi-infinite medium, which is very different from these superluminal phenomena caused by reshaping¹⁶ or tunneling¹⁸ effects in a finite region. Our results confirm that the wave energy propagation velocity in a real system cannot exceeds the speed of light in vacuum, which is a restriction derived from special relativity.