Characteristic wave velocities in spherical electromagnetic cloaks

We investigate the characteristic wave velocities in spherical electromagnetic cloaks, namely, phase, ray, group and energy-transport velocities. After deriving explicit expressions for the phase and ray velocities (the latter defined as the phase velocity along the direction of the Poynting vector), special attention is given to the determination of group and energy-transport velocities, because a cursory application of conventional formulae for local group and energy-transport velocities can lead to a discrepancy between these velocities if the permittivity and permeability dyadics are not equal over a frequency range about the center frequency. In contrast, a general theorem can be proven from Maxwell's equations that the local group and energy-transport velocities are equal in linear, lossless, frequency dispersive, source-free bianisotropic material. This apparent paradox is explained by showing that the local fields of the spherical cloak uncouple into an E wave and an H wave, each with its own group and energy-transport velocities, and that the group and energy-transport velocities of either the E wave or the H wave are equal and thus satisfy the general theorem.


Introduction
Invisibility cloaking has received considerable attention since the theory was first presented in [1] and experimentally demonstrated to within an approximation in [2]. The approach presented in [1] provides a conceptually simple method for designing a cloak. One can imagine that the trajectories of electromagnetic waves passing through a region of warped space must conform to the local metric. Once the desired trajectory is determined through a transformation applied to Cartesian straight trajectories, the differential operators in Maxwell's equations in the transformed space lead to space-dependent coefficients that can be reinterpreted in terms of constitutive relations of an anisotropic, inhomogeneous medium. In [3], it was shown that cloaking can be reformulated as a boundary value problem with a single first-order Maxwell differential equation for linear anisotropic media, leading to explicit formulae for the relative permittivity and permeability dyadics and fields of ideal spherical and circular cylindrical annular cloaks.
Causality-energy conditions imply that electromagnetic incident fields with a finite bandwidth cannot be perfectly cloaked [3]. Consequently, it is desirable to evaluate the performance degradation of electromagnetic cloaks as a function of the frequency for permittivity and permeability dyadics that satisfy realistic causality-energy relations [4]- [8]. Toward this end, we investigate some characteristic wave velocities of spherical electromagnetic cloaks as determined by the Poynting vector, the local differential propagation constant, and the phase velocity. Of particular interest are the group and energy-transport velocities, which provide information on the potential cloaking bandwidth.
Explicit formulae are provided for the phase velocity and the Poynting vector inside an ideal spherical cloak illuminated by a plane wave. It is shown that the 'complex' Poynting vector is real-valued and divergence-free inside the cloak, and that it is aligned with the 'compressed' rays obtained with the radial transformation function [9]. On the other hand, the phase velocity is not aligned with the Poynting vector, thus leading to the definition of a 'ray velocity' as the phase velocity along the direction of the Poynting vector.
The calculation of the group and energy-transport velocities requires a realistic frequency dependence for the permittivity and permeability dyadics. Since the frequency dependences of 3 these two dyadics may be significantly different in an anisotropic material, it is useful to assume physics-based forms for these functions in order to investigate the degradation in cloaking as a function of frequency about a center frequency ω 0 at which the relative permittivity and permeability dyadics are equal. We calculate the local group velocity v g = ∇ k ω(k) and the energy-transport velocity defined as the time-average Poynting vector S(k) divided by the energy density U (k) [v e = S(k)/U (k)] for a general spherical cloak. If the variation with frequency of the relative permeability and permittivity dyadics are the same, we find that the group and energy-transport velocities are identical (v g = v e ). However, if the variation with frequency of the relative permeability and permittivity dyadics are different, it appears that the group and energy-transport velocities are not necessarily equal even at the center frequency where the relative permeability and permittivity are the same [4].
This apparent discrepancy between the values of the local group velocity and the energytransport velocity is paradoxical in view of the general theorem provable from Maxwell's equations for linear, lossless, frequency dispersive and anisotropic material, namely that the local group velocity and the energy-transport velocity are identically equal. In order to resolve this apparent paradox, we provide an especially transparent derivation of this general theorem and examine its application to the material of dispersive spherical cloaks [10].
The paper is organized as follows: section 2 summarizes the results for the phase and ray velocities in spherical cloaks at their operational frequencies. Section 3 provides an especially transparent proof of the equality between the local group velocity and the energy-transport velocity in linear, lossless, frequency dispersive and bianisotropic material. In section 4, the electromagnetic fields of the spherical cloaks are shown to uncouple into E and H waves, and the local group and energy-transport velocities are proven equal for the individual E and H waves. Numerical results for the ray, group and energy-transport velocities are presented in section 5 and the paper ends with a concluding section 6.

Ray and phase velocities in spherical cloaks
The boundary value formulation of electromagnetic cloaking presented in [3] is based on the requirements that the cloaking occurs for all possible incident fields, the cloaks have continuous tangential components of E and H fields across their outer surfaces, and the normal components of the D and B fields are zero at the inner material surfaces of the cloaks. The tangentialfield boundary conditions at the outer surface of a cloak ensure zero scattered fields, and the normal-field boundary conditions at the inner surface of a cloak are compatible with zero total fields inside the interior cavity of the cloak. For a cloak consisting of a spherical annulus of anisotropic material with inner radius a and outer radius b, these boundary conditions lead to the following expressions for the relative permittivity and permeability dyadics at the operational frequency ω 0 where 0 and µ 0 are the free-space permittivity and permeability, respectively, r is the radial coordinate of a spherical coordinate system (r, θ, φ) with origin at the center of the spherical cloak, f (r ) is the 'compressed radial coordinate function' satisfying the boundary conditions f (a) = 0 and f (b) = b, and the prime superscripts denote differentiation with respect to r .

4
The field distributions associated with the relative permittivity-permeability dyadic in (1) are where (E inc , H inc ) are the incident fields. The particularly simple choice for the radial function yields the spherical cloak of Pendry et al [1].

Poynting vector and phase velocity
It is interesting to investigate the simple case of a plane wave impinging on the spherical cloak. It is not restrictive to assume that the plane wave travels in the −ẑ-direction with electric field polarized alongŷ, that is where ζ 0 is the free-space impedance and an exp(−iω 0 t), ω 0 > 0 time dependence has been assumed and suppressed. Inserting these equations into (2), one obtains From these fields, a simple expression can be derived for the time-average Poynting vector as The 'complex' Poynting vector, E × H * /2, inside the cloak is real, and this is true whenever the incident field exhibits a local plane-wave structure in the region captured by the space compression. Furthermore, it is easily verified that ∇ · S = 0 everywhere, a result that confirms that there is no power loss at any point of the cloak. The unit vectorŝ in the direction of the Poynting vector within the material of the spherical cloak is given bŷ This unit vectorŝ defines the tangent to the power-flow rays at each point in the material of the cloak and can be interpreted in terms of a generalized ray-direction, thus leading to pictures that resemble ray-path propagation in inhomogeneous isotropic media. However, conventional ray interpretation is not fully possible because of the anisotropy of the medium. The generalized ray-paths inside the cloak are plotted in figure 1(a) from (7) for the case b = 2a. It is apparent that the rays penetrate the cloaking material while avoiding the central free-space cavity.
Within the transformational optics framework (see, for instance [1]), the power-flow rays are directly defined from a mapping of the free-space incident ray-fields to the incident ray-fields 'compressed' within the spherical cloak material by the coordinate transformation function f (r ). We shall now show that the rays obtained from this mapping are identical to the powerflow rays obtained from the Poynting vector. Indeed, the equation of the incident ray-field at a distance ρ 0 from the z-axis is r = ρ 0 /sinθ , while the equation of the 'compressed' ray-field within the material of the spherical cloak is f (r ) = ρ 0 /sin θ from which we obtain The unit tangent vectort at a generic point in the cloak iŝ t =θ r dθ +rdr Comparing (9b) with (7) shows thatt =ŝ, that is, the unit tangent vector at any point of the rays formed by simply mapping the original rays through the field transformation function f (r ) is equal to the unit vector in the direction of the Poynting vector at that point. We have proven this for each straight-line ray of an incident plane wave illuminating the spherical cloak. However, because the incident rays of any source (lying outside the cloak) illuminating the spherical cloak are straight lines in free space, each incident ray can be considered locally as part of a plane wave with an unlimitedly high frequency. Thus, the proof holds for any incident ray illuminating the spherical cloak. In other words, the rays defined by the Poynting vector in a spherical cloak illuminated by an arbitrary incident field are identical to the rays formed by mapping the incident rays through the compression function f (r ).
The equiphase wavefronts are obtained by equating the argument of the exponential in equations (5) to a constant: f (r )cos θ = constant, so that f (r )cos θ dr − f (r ) sin θ dθ = 0. The unit tangent vector to these equiphase wavefronts is then τ =θ r dθ +rdr which is plotted in figure 1(b). Comparing (10) with (7) shows that the tangents to the wavefronts are rotated 90 • with respect to the Poynting vector lines (power-flow rays). This occurs at any plane of constant φ and holds for any compression function f (r ).
The local differential propagation constant normal to the equiphase wavefront surfaces is given by and its direction is shown in figure 1(c). It is worth noting that the local normal to the wavefronts is not aligned with the Poynting vector, in contrast to what happens for conventional rays. The magnitude of the local differential phase velocity normal to the wavefronts is where c is the free-space speed of light. The magnitude of the phase velocity along the Poyntingvector rays is therefore Along the ray path inside the cloak, the generalized rays propagate with an average phase velocity greater than the speed of light and recover, as they exit, the same phase as the plane wave along the straight path. The latter statement can be easily verified by considering an incident ray at a distance ρ 0 from the z-axis, as shown in figure 2. In the absence of the cloak, the time delay accumulated by the ray along the path from the point A to the point B is where 2z 0 is the distance between A and B. On the other hand, in the presence of the cloak, the time delay is Hence, as expected, all the rays travel from one side of the cloak to the other in the same time as they would take in the absence of the cloak.

Group and energy-transport velocities in bianisotropic media
Consider Maxwell's equations in k-ω space for lossless, homogeneous, spatially nondispersive, source-free media with bianisotropic constitutive relations in which the exp[i(k · r − ωt)] space-time dependence has been suppressed, where the rectangular components of the propagation vector k can take any values on the real line, (−∞, −∞, −∞) < (k x , k y , k z ) < (+∞, +∞, +∞). In a lossless medium, the permittivity, permeability and magnetoelectric constitutive dyadics obey the relations Elimination of all the fields but E from (16) and (17) from which the dispersion equation follows: where k is the antisymmetric dyadic defined as Because the medium is lossless and k is a real propagation vector, energy conservation requires that the solutions ω(k) to (20) be real. Also, it can be proven [12] that for lossless, homogeneous, spatially non-dispersive, bianisotropic material that is also reciprocal, solutions for a given ω come in pairs with propagation vectors k and −k * , which equals −k for the real k propagation vectors considered here in three-dimensional (3D) Fourier k-space.
Depending on the values of the constitutive dyadics, more than one characteristic value of each positive or negative ω(k) may satisfy (20) and each of these characteristic values generally has a different associated characteristic field solution to (16) and (17). If this is the case, that is, there exists more than one characteristic wave solution to (16) and (20), then in the following derivations, ω(k) refers to any one of these characteristic values and the fields are those associated with that particular characteristic value.

Wavepacket and group velocity
The real space-time fields in the medium can be expressed as a four-fold integral consisting of a 3D Fourier transform of the fields in k space and an analytic Fourier transform over positive frequencies ω [13, chapter 5]. Moreover, for any one characteristic value of ω(k) > 0 obtained from the dispersion equation (20), the frequency Fourier transform reduces to this one discrete value and we have for the electric field with dk = dk x dk y dk z and the integration limits k ± are (−∞, −∞, −∞) < (k x , k y , k z ) < (+∞, +∞, +∞). A wavepacket is defined by assuming that the wave propagation vectors are concentrated about a central propagation vector k 0 such that E(k) ≈ 0 unless |k − k 0 |/|k 0 | = | k|/|k 0 | 1. Assuming ω(k) is expandable about k 0 , we have so that E(r, t) in (22) can be approximated as or simply where ω 0 = ω(k 0 ) and g(r) is an envelope function that varies slowly over a distance equal to the wavelength λ 0 = 2π/|k 0 | in the medium. The velocity of this envelope or wavepacket is called the group velocity and is determined from (24b) as which is not necessarily in the same direction as the phase velocity v p = (ω 0 /|k 0 |)k 0 . As previously explained, if the material is reciprocal, the group velocities in the ±k 0 directions are equal and thus only values of (k 0x , k 0y , k 0z ) > 0 as well as ω(k) > 0 need to be considered for reciprocal material.

Proof that group and energy-transport velocities are equal
Taking the differential of (16a) and (16b) dotted into H * and E * , respectively, one obtains By subtracting these two equations and then using the constitutive and lossless relations in (17) and (18) for bianisotropic media leads to where is the time average Poynting vector (power flow per unit area) and is the energy density in a lossless, homogeneous, frequency dispersive (but spatially nondispersive), bianisotropic medium [14]. Since (27) holds for all dk, it follows that where v e (k) = S(k)/U (k) is the 'energy-transport velocity' for a wavepacket with phase velocity in the k-direction. The relationships in (25), (30) and (31) can be summarized in the one extended equation which says that the group and energy-transport velocities of a wavepacket are equal in a lossless, homogeneous, frequency dispersive (but spatially non-dispersive), bianisotropic medium.

Evaluation of group and energy-transport velocities in a spherical cloak
In order to evaluate the group and energy-transport velocities in a spherical cloak, specific frequency dependencies must be assumed for the relative constitutive dyadics. These frequency dependencies may be significantly different for the permittivity and permeability dyadics, but they are assumed to be equal in accordance with (1) at a center frequency ω 0 . Reasonable functional forms for these frequency dependencies will be explicitly given in the next section.
Here, we simply let the relative permittivity and permeability dyadics have the general frequency dependence where (ω 0 )/ 0 = µ(ω 0 )/µ 0 such that but they are not necessarily equal for ω = ω 0 . Note that r , s , µ r , µ s are elements of the relative permittivity and permeability dyadics. The expressions in (25) and (31) for the group velocity and energy-transport velocity have been derived under the assumption that the medium is homogeneous. In order to apply these equations to the inhomogeneous material of the spherical cloak, we let the local values of (ω)/ 0 and µ(ω)/µ 0 in (33) form part of an infinite homogeneous medium, insert those values into (20), and solve for ω(k). (The group and energy-transport velocities evaluated under this assumption of 'local homogeneity' give the 'ray-optics approximation' to the velocity of an actual pulse propagating in an inhomogeneous spherical cloak-an approximation that becomes more accurate with less variation per wavelength in (ω)/ 0 and µ(ω)/µ 0 [15, sections 1.6 and 1.7]. The issue that we are addressing here in section 4 is the discrepancy that occurs between the group and energy-transport velocities if it is not realized that the fields in the cloak separate into uncoupled E and H waves. This discrepancy exists independently of the accuracy of the ray-optics approximation of local homogeneity that allows us to apply (25) and (31).) Letting the local spherical components (k r , k θ , k φ = 0) coincide with the rectangular components of k of the wavepacket in the associated infinite homogeneous medium, we find that the dispersion equation (20) can be locally factored as where The individual dispersion equations F E (k, ω) = 0 and F H (k, ω) = 0 admit two, generally different uncoupled solutions, ω = ω E (k) and ω = ω H (k), respectively, corresponding to a local E wave (H φ = 0) and a local H wave (E φ = 0); specifically The individual group velocities of the E and H waves in (37) and (38) can be determined from the identities Unless (41) is satisfied, the group and energy-transport velocities in (42b) and (43) are not equal and neither of them give the group-energy velocity in (40) of the E wave or the H wave. Failure to recognize that the spherical-cloak material supports two uncoupled wavepackets and that these wavepackets can have different group-energy velocities if the variation of µ(ω)/µ 0 and (ω)/ 0 do not satisfy (41) for ω = ω 0 , even though µ(ω 0 )/µ 0 = (ω 0 )/ 0 , leads to the imperfect results in (42b) and (43) for the group and energy-transport velocities in the spherical cloak, rather than the correct results in (40).

Numerical results
In order to numerically evaluate the group and energy-transport velocities inside the spherical cloak, we shall assume the same realistic frequency variations for the relative permittivity and permeability dyadics as those adopted in [4]: • s (ω) and µ s (ω) can be assumed to be slowly varying with frequency since s (ω 0 ) = µ s (ω 0 ) > 1, and thus approximately constant • a Drude model is assumed for r (ω), namely • a Lorentz model is assumed for µ r (ω), namely µ r (ω) = 1 + Note that these values of r , s , µ r , µ s satisfy (34) at ω = ω 0 . The expressions in (44)-(46) for the constituent parameters have been used together with (3) and (40) to numerically calculate the values of the group-energy velocities at the frequency ω 0 inside a spherical cloak with b = 2a. Figure 3 shows the values obtained for  different values of the incident ray distance ρ 0 from the z-axis. It is apparent that the groupenergy velocities for the E and H waves have similar behaviors but different values inside the cloak. The corresponding values of the ray velocity given in (13) are also depicted on the same plot. Finally, bidimensional plots of the ray velocity and the group-energy velocities for E and H waves are shown in figures 4(a)-(c), respectively. The generalized ray trajectories [ŝ(k 0 )], which are the same for both the E and H waves, are also shown in these figures.
As the plots in figures 3 and 4 indicate, there is an infinitely large pulse time delay for the ray that is directed toward the center of the cloak, a result found previously by Chen and Chan [4].

Conclusion
Characteristic phase, ray, group and energy-transport velocities have been determined in spherical electromagnetic cloaks with material constitutive parameters having realistic causal frequency variations. After providing explicit expressions for the phase and ray velocities, it is confirmed that the ray travel time across the space occupied by the cloak is the same with and without the cloaking material. Particular attention is given to the group and energytransport velocities. The general theorem stating that the group velocity [∇ k ω(k 0 )] equals the energy-transport velocity [S(k 0 )/U (k 0 )] in linear, lossless, homogeneous, frequency dispersive, source-free, bianisotropic material applies to a single characteristic local wavepacket in that material. If the dispersion equation allows for more than one characteristic solution and thus more than one characteristic positive frequency ω(k) and wavepacket, the general theorem applies to each separate characteristic frequency and wavepacket but not to the combined dispersion equation and the combined fields of the characteristic waves. The discrepancy between local group and energy-transport velocities in the anisotropic material of spherical cloaks, that arises from a cursory application of the formulae for these velocities, is traced to the failure to recognize that the wavepacket solution in the spherical-cloak material uncouples into E and H characteristic waves with different characteristic frequencies if the frequency variation of the relative permittivity and permeability dyadics are not the same. A necessary