Envelope Dyadic Green’s Function for Uniaxial Metamaterials

We introduce the concept of the envelope dyadic Green’s function (EDGF) and present a formalism to study the propagation of electromagnetic fields with slowly varying amplitude (EMFSVA) in dispersive anisotropic media with two dyadic constitutive parameters: the dielectric permittivity and the magnetic permeability. We find the matrix elements of the EDGFs by applying the formalism for uniaxial anisotropic metamaterials. We present the relations for the velocity of the EMFSVA envelopes which agree with the known definition of the group velocity in dispersive media. We consider examples of propagation of the EMFSVA passing through active and passive media with the Lorentz and the Drude type dispersions, demonstrating beam focusing in hyperbolic media and superluminal propagation in media with inverted population. The results of this paper are applicable to the propagation of modulated electromagnetic fields and slowly varying amplitude fluctuations of such fields through frequency dispersive and dissipative (or active) anisotropic metamaterials. The developed approach can be also used for the analysis of metamaterial-based waveguides, filters, and delay lines.

function methods or the full-wave numerical simulation methods, is in a great reduction of the computation time, because the analytical expressions for the EDGF matrix elements are derived as explicit functions of the space and time coordinates. The uniaxial MMs are also known for advantageous optical properties for sensing 28 , nonlinear optical applications 29 , and spontaneous emission control 30 . In particular, in extremely anisotropic uniaxial MMs the beam diffraction effects can be suppressed, which creates new opportunities for the flexible control of the RHT in MMs.
Note that the above-mentioned effects depend on the relation between the phase and group velocities in a MM. In passive media, the sign of the phase velocity, ω = v k / ph , can be chosen with respect to the direction of the group velocity, ω = ∂ ∂ v k / gr , which is typically selected to be always positive. This is rooted in the fact that, in low dispersive passive media, the directions of the energy and information transfer and the group velocity coincide. However, in active media, the signal (i.e., information) propagation direction does not necessarily coincide with the apparent energy flow direction, because such media can produce responses more energetic than the incoming signals. On the other hand, the information propagation direction in any media is fixed by causality, which means that the sign of the group velocity can be uniquely determined with respect to this direction. Under this convention -which we adopt in our work -the sign of the group velocity can vary depending on the material properties. Moreover, in dispersive media with loss, the group velocity can become complex. Although the real part of the complex group velocity can be still related to the envelope's propagation speed, the interpretation of the imaginary part of the group velocity varies among different authors 31 .
In some situations, the unusual dispersive properties of the MMs result in superluminal or subluminal group velocities, in which case this quantity can be higher 32 or extremely lower than the speed of light in a vacuum 33 . It can also become zero or negative [34][35][36] . These effects in MMs, as well as the negative refractive index, are well-accommodated within the framework of the causality principle [37][38][39][40] . In particular, the superluminal effects are due to interference between the nearby frequency components in a narrow-band signal that has experienced propagation through a medium with anomalous dispersion 41 . Such effects are essentially kinematic and are not associated with any superluminal exchange of information and (or) energy between a signal source and a receiver. As is known 37 , the transfer of information is associated with the propagation of the envelope's precursor, rather than with the envelope's maximum. The former always happens at a speed ≤ v c. In recent decades, superluminal and subluminal group velocities have attracted attention in nonlinear optics 42 , quantum communication 43 , photon controlling and storage [44][45][46] , precision sensing 47 , high-speed optical switching 48,49 , broadband electromagnetic devices and delay compensation circuits in ultra-high-speed communication systems 50 , and in high resolution spectrometers 51 . In this work, we examine the EDGF formalism on the hyperbolic MMs and media with gain, in order to demonstrate the applicability of the method to the wave processes in such media, including the negative refraction and focusing and the exotic superluminal processes in active media.
This article is organized as follows. We start from introducing the slowly varying complex amplitudes (SVA) of the electric and magnetic fields and writing the macroscopic Maxwell equations and the material relations in terms of these quantities. The EDGF is then found as the solution of the resulting 6-dyadic dynamic equation with a band-limited point source. Next, we obtain the EDGF matrix elements for the uniaxial anisotropic MMs and develop the paraxial approximation for the EDGF. In numerical examples, we apply the EDGF method to study characteristic effects associated with the propagation of spatio-temporally modulated signals in MMs (negative refraction) and in active media (superluminality), and confirm the applicability of the EDGF approach to these problems. We conclude with discussion of the obtained results and outline possible generalizations of the developed formalism.

Results
envelope dyadic Green's function technique. In order to investigate the propagation of the EMFSVA in a dispersive medium, we consider the Maxwell equations for the time-dependent electromagnetic fields and sources written as follows , the electric (magnetic) charge and current densities, respectively, are related by the continuity equation, . The electromagnetic constitutive equation is given by¯ζ where the dyadic integro-differential operator ¯ζ ε = (μ¯) represents the dispersive anisotropic permittivity (permeability) of the material, which relates Ξ = D B ( ) to = f E H ( ). The time-dependent electromagnetic field can be expanded into the monochromatic spectral components as follows On the other hand, considering the time-harmonic fields with SVAs, it is acceptable that the most spectral energy is concentrated in a narrow band Δω around ω 0 , the carrier frequency, with ω ω Δ  0 . Thus, we can define the EMFSVA, t f ( ) m , as follows where c.c. is the abbreviation for the complex conjugate of the previous term and ω ω Ω = − 0 . We also assume that ζ ω , the Fourier transform of ζ , has appreciably smooth behavior in the narrow band Δω, which is a reasonable assumption for applications of our interest. So, expanding ζ ω around ω 0 and keeping only the two first terms, we obtain 0 0¯w and | ω 0 denotes that the related quantity is computed at ω 0 . Using Eqs. (5-7), we obtain Ξ ∂ t as follows , and we have to ignore the second and higher orders of ∂ t in expressions involving this operator. With having these tools at hand, we proceed to the Green's function technique to solve the Maxwell equations written in terms of the EMFSVA.
Starting from Eqs. (1) and (2) for the time-dependent electromagnetic fields and sources and using Eqs. (3)(4)(5)(6)(7)(8)(9), the Maxwell 6-vector equations for the EMFSVA assume the following operator form: In the above equation, F m , the SVA of the 6-vector field, J m , the 6-vector current density, and Ô , the electromagnetic operator, are given by where Ī is the identity dyadic. In Eq. (11), the differential operators ī ( ) ith the second order derivative term ignored, in line with the assumption of the SVA. We may introduce the envelope 6-dyadic Green's function (EDGF), t t G r r ( , ; , ) ′ ′ , for Eq. (10) such that where I is the identity 6-dyadic, and δ ⋅  ( ), similar in role to the Dirac delta δ ⋅ ( ), emphasizes that we consider the slowly varying sources and fields. Respectively, δ  t ( ) has a finite duration and δ | | < ∞  t ( ) . Equation (13) means that the EDGF is obtained by inverting the operator Ô .
It is more convenient to work in the Fourier domain when the material parameters are uniform in both space and time. In this case, t t t t G r r Gr r ( , ; , ) ( , ) ′ ′ ≡ − ′ − ′ , and we can define g k ( , ) Ω , the Fourier transform of the EDGF, by the following relation Recalling Eq. (13) and using that, in the Fourier space, ∂ = − Ω i t and ∇ = ik, where k is the wave vector, we find that the 6-dyadic equation for g satisfieŝˆ¯¯¯ where we have used operator-like notation for Ω D such that . Knowing the dispersive constitutive parameters of the medium and inverting the operator Ô in Eq. (15), we can obtain g and, in turn, the EDGF. Afterward, for any given forms of the 6-vector source functions J m with the frequency spectra fitting the interval ; 2 2  In the next section, we shall find the matrix elements of the EDGF for uniaxial anisotropic media. eDGf matrix elements. Here, we apply the EDGF technique for reciprocal uniaxial media to obtain the matrix elements. For such media we can write xx yy tˆˆˆˆ and ˆÎ zz z = being the projection dyadics in the transverse and axial (the main axis) directions, respectively, along which the electromagnetic responses of the medium are different. In order to obtain the EDGF of the uniaxial medium, we impose this property of constitutive parameters to decompose the operator Ô in Eq. (15) as (19) where 0 is the null dyadic. When writing the matrix elements in Eq. (19), we keep in mind that Next, it is more convenient if we replace the second row with the third one and also the second column with the third one in the matrix representation of Ô in Eq. (19) so that the tangential and axial components of the constitutive parameters of the medium take place in separate blocks. When doing this, it is also necessary to respectively rearrange the dyadic components of g and I in Eq. (15). After doing this and taking into account that I I k k t t z t × = ×¯, the EDGF is expressed through the inverse of the reordered matrix of Ô as follows (20) Representation (20) allows for a straightforward splitting of the fields into a pair of orthogonal polarizations: The transverse-electric (TE or s-) polarization with vanishing axial component of the electric field, and the transverse-magnetic (TM or p-) polarization with vanishing axial component of the magnetic field. In order to perform such splitting, Ī t in Eq. (20) is expanded as ¯= Then, after a rather tedious but straightforward dyadic algebra, the matrix in Eq. (20) can be inverted and the following result obtained:  www.nature.com/scientificreports www.nature.com/scientificreports/ where the 12 independent non-vanishing components of g are expressed as follows (here and thereafter we use shorthand notation α β [ ] for respective selection of either α or β, where α and β are isolated terms or groups of indices): where, from Eq. (16), where the index l is either t or z, and As can be seen from Eq. (24), there are two poles in k z for each polarization: p κ ± for the TM case and, similarly, s κ ± for the TE one. Physically, the poles κ = + k z ps [ ] correspond to the waves propagating in the halfspace − ′ > z z 0, and the poles with κ = − k z ps [ ] correspond to the waves propagating at − ′ < z z 0. Based on Eq. (24), we can sort out the components of g into two groups which correspond to the waves of pand s-polarization where N p s [ ] are formed by the corresponding terms in the numerators of Eq. (24). Respectively, when taking the inverse Fourier transform of g as given by Eq. (14) and considering the integral over dk z , we get two categories of residues 1 is the sign of Z. Next, before taking the integral over Ω d , we recall the approximation of the SVA, and expand γ ± A , 1 and κ γ around the point Ω = 0 and ignore Ω O( ) 2 terms so that  for the p-[s-] polarization. It should be noted that ignoring the second order terms in Eqs. (29) and (30) is reasonable for the frequency intervals of our interest, where the group velocity dispersion effects are relatively weak. In Methods, we discuss this with more detail and also present a closer look at the group velocity for a propagating envelope in an active and dispersive medium.
With these approximations at hand, when performing the integration over Ω we obtain denotes the spherical Bessel function of the first kind and the n-th order. Therefore, from Eq. (14), we obtain the EDGF, τ G R ( , ), in the following form , 1 , 1 0 , 1 1 As can be inferred from Eqs. (34) and (35), besides the usual effect of the imaginary part of the propagation factor, κ γ Im( ) 0 , the imaginary part of the complex group velocity, γ V Im( ) g , is responsible for an additional change in the envelope's magnitude during propagation.
In some special cases, the remaining integration over dk t in Eq. (34) can be performed analytically. In the following subsection, we consider such a special case of paraxial propagation. The representation (34) is most suitable for the calculation of fields of sources with a known spatial spectrum in the transverse plane, e.g. in near-field radiative heat transfer problems. Considering arbitrary source vectors and using Eqs. (17) and (34), we can investigate the propagation of the EMFSVA in any unbounded anisotropic dispersive media.
The present formalism can be extended to contain interface effects, which will make it applicable to investigate the propagation of the EMFSVA via multilayer media. In what follows, we consider a special case of such media when the neighboring layers are (approximately) impedance matched. eDGf for paraxial propagation. Let us consider the case when the envelope propagation happens dominantly along the anisotropy axis. This case is typical for extremely anisotropic uniaxial MM in which ε ε | | | |  z t and (or) μ μ | | | |  z t . Indeed, the propagation factor κ γ 0 from Eq. (31) can be expressed as On the other hand, from Eq. (24) it is seen that in this case the g tt component dominates over the g tz , g zt and g zz components, so that the EDGF is dominantly transverse which implies that the wave energy propagates dominantly along the z-axis.
Therefore, when calculating the integral over dk t in Eq. (34), we will not make a big mistake if we evaluate the γ ± T , 1 term of Eq. (34) at → k 0 t while using the expansion (37) in the exponential terms. In this approximation, the integration over dk t can be performed analytically, which results in the following expression for the tt-block of the paraxial EDGF: The integral over ξ d in Eq. (38) can be taken in a closed form after the dyadic differential operators t γ± have acted on the exponential term. The integral representation (38) is especially handy when taking the convolution integral in Eq. (17) with the source currents having Gaussian profiles in the xy-plane. . Under this condition, practically all source spectral power is concentrated within the frequency interval of width ω Δ . Thus, the SVA of the 6-vector source current density reads Such a source creates the electromagnetic field in the 0-th layer which propagates in both > ′ z z and < ′ z z directions. When this field reaches the interface = z z 1 between the 0th and the next layer, it excites the fields in the next layer, and so on. If the characteristic wave impedances of the neighboring layers are mismatched, the reflected field will also appear, which can propagate to the other interface, be partially reflected again, etc.
We reserve the study of such multiple reflections in the EDGF context for a future work. Here, we assume that the neighboring layers are approximately impedance matched at the frequencies close to ω 0 , and the reflections may be neglected. Note that this does not mean that the layers must be made of the same materials, or that the materials must have the same dispersion. The impedance match condition for the case of the paraxial propagation considered in this section, means that the material parameters of the l-th and + l ( 1)-th layer satisfy Under this assumption, the EMFSVA created by the source (42) in the 0-th layer can be obtained from Eqs. (17) and (38) and expressed in the following form, after evaluating all involved integrals:¯¯¯    , etc. are expressed through the components of the material dyadics of the 0-th layer, ε 0 and μ 0 .
In order to find the fields in the next layer located at ∈ z z z ( , ) 1 2 (the 1st layer), we introduce an equivalent Huygens source (a pair of electric and magnetic surface currents) placed at = z z 1 (where > ′ z z 1 ) defined as The field in the layer ∈ z z z ( , ) 1  numerical examples. Negative refraction and focusing by uniaxial MM with hyperbolic dispersion. When the transverse and longitudinal components of the permittivity dyadic of a MM have opposite signs, the isofrequency curves for the p-polarized waves are hyperbolas. It is known that the p-polarized light refracts negatively when impinging at the interface of a conventional material and such a hyperbolic MM 53,54 . In the following numerical example, we study the implications of this phenomenon on the paraxial propagation of the Gaussian envelopes considered earlier. We shall confirm that the EDGF formalism correctly predicts focusing of a diverging p-polarized Gaussian beam by the hyperbolic MM. At near-infrared frequencies, a hyperbolic MM can be realized, e.g. by embedding vertically aligned metallic nanowires into an isotropic dielectric host. For the sake of a numerical example, here we consider a MM formed by golden nanowires embedded into alumina substrate. By using the Maxwell-Garnett effective medium theory (EMT) for a uniaxial MM formed by such nanowires, the following expressions for the effective transverse and axial permittivities can be obtained 55 : where ε h and ε m are the dielectric permittivities of the host material (Al 2 O 3 56 ) and the plasmonic metal (Au), respectively, and f is the nanowires volume fraction. The relative permittivity of gold at near-infrared frequencies follows the Drude dispersion model 57 where ω = . [Eq. (49)]. The results of the paraxial EDGF-based calculations for this structure for the p-and s-polarized source currents are displayed in Fig. 1. As one can see, the initially diverging p -polarized Gaussian beam, after refracting negatively at the interface λ = z/ 10 0 , is focused at λ = z/ 20 0 , the middle point of the hyperbolic MM layer, and then it diverges again. After reaching the second interface at λ = z/ 30 0 , the beam undergoes another negative refraction and is focused inside the second dielectric layer at the point λ = z/ 40 0 . On the contrary, the s-polarized beam does not experience any refraction at the MM interfaces and simply diverges. Note that the scales on the x-axis and the z-axis in Fig. 1 are different, so that the field profile along z is compressed in comparison with that along x.
In order to explain how our EDGF formalism is able to reproduce these phenomena in the paraxial approximation, let us consider the expression for the square of the effective beam width, σ . This source creates the EMFSVA propagating through the three media in the > z 0 direction. Because the relative permittivities and permeabilities of the layers are rather close to unity, the layers are well impedance matched and we can apply the theory developed earlier.
The relative permittivity of the 132 Xe gas with inverted population follows the Lorentzian dispersion with negative oscillator strength 34 : , the maximum of the original pulse passes forward and the field acquires a noticeable value at the interface of the first and the second layers and, at the same time, we can see that the amplified field in the active layer forms a sharp peak and penetrates into the third layer. We can see how the back propagating pulse is formed in the gain medium and how it interacts with the primary pulse which propagates forward in the first layer at the same interface at ω = t 260, 300, p and 360. Finally, at ω = t 500 p , we can see that the main pulse has left the two layers, however, a small effect of its back edge is still present at the interface of the first and the second layers. It should be noted that at ω = t 260 p the maximum of the pulse exits the second layer earlier than it would do if it had traveled through an equal distance of air, i.e. it appears superluminal. This happens due to the action of the gain medium on the electromagnetic field with a Gaussian-shaped envelope which lacks a definite turn-on moment.

Discussion
In this work, we have presented the EDGF concept and the related formalism which are applicable for studying the propagation of the amplitude fluctuations of quasi-monochromatic electromagnetic field in anisotropic dispersive media. The developed formalism is aimed to be used in future works to model the dynamics of RHT in such media, in particular, in uniaxial MMs 22 , however, it is equally applicable to the analysis of narrow-band signal propagation in these MMs. Although not considered in this article, the present formulation allows also for a direct finite-difference time-domain-based numerical solution of the dynamic equations for the EMFSVA, which can be useful, e.g. when studying propagation of wave packets in layered MMs.
Note that, in contrast to the conventional dyadic Green's functions (DGF), our method incorporates the SVA approximations into the linear operator and the source of the Green's equation. Thus, the EDGF is a wavelet-type dyadic function that describes excitation and propagation of the band-limited amplitude and phase fluctuations of the quasi-monochromatic electromagnetic fields. These fluctuations propagate with a velocity that has the meaning of the group velocity. This makes the EDGF approach distinct from the DGF formulations that deal with delta-functional sources in both space and time and result in signal fronts propagating with the speed of light in www.nature.com/scientificreports www.nature.com/scientificreports/ vacuum, ε μ = c 1/ 0 0 , due to the fact that in any physical material ¯( ) Moreover, it can be shown that the EDGF approach inherits from the SVA approach an ability to deal with materials with slow and weak nonlinearity, in which the constitutive parameters depend on the amount of electromagnetic energy passed through the medium 58 .
There is a significant body of literature on the standard DGF in anisotropic and bianisotropic media [59][60][61][62][63][64][65] . While in the most general case of bianisotropic medium there is no closed-form representation for the DGF, in the case of uniaxial magnetodielectrics, the time-harmonic DGF can be expressed through a pair of scalar electric and magnetic Green's functions. When interested in the EMFSVA in such media, a possible approach could be to start from such closed-form representations and expand the frequency-dependent parameters (such as the material dyadics, the propagation factors, etc.) in these relations into series around the carrier frequency ω 0 , with ω ω Ω = − 0 being the small parameter.
Although this would constitute a tractable approach for some cases of uniaxial and bianisotropic media, in this article we have developed a conceptually different method, which can be later extended to media with more general material relations, including nonlinear media (e.g.by using Green-Volterra series 66 ). By incorporating the SVA approximations into the material relations and the Maxwell equations, we have formulated the 6-dyadic equation for the EDGF in anisotropic dispersive media and obtained the matrix elements of the EDGF for the uniaxial magnetodielectric media in the Fourier and the configuration spaces. In the case of paraxial propagation, the EDGF for uniaxial media can be written in a closed form, resulting in a formulation similar to the Gaussian beam-based paraxial approximation in optics.
In order to test the ranges of applicability of our formalism, we have considered the propagation of the EMFSVA through a hyperbolic MM and an active medium, sandwiched between two passive dielectric layers. In the former case, the paraxial propagation of the EMFSVA leads to the well-known negative refraction effect 53,54 , while in the latter case one can observe the peculiar phenomenon of negative group velocity associated with superluminal propagation 34 . We have modeled the propagation of the EMFSVA through such media by employing the effective Huygens sources at the interfaces of the neighboring layers. It is worth to mention that the paraxial propagation of the EMFSVA can be used to study the near-field thermal radiation effects in extremely anisotropic media, e.g. in arrays of aligned carbon nanotubes 67 .
We have demonstrated with numerical examples that the developed formalism correctly models negative refraction and focusing by hyperbolic MMs and is also applicable to exotic effects in optical media with inverted population, such as the negative group velocity and superluminality. The group velocity obtained in our formalism agrees with the standard definition known from the literature and results in similar behaviors. The considered examples confirm the applicability and validity of the EDGF approach developed in this paper. for the considered example. This value can be considered sufficiently small. On the other hand, in the regions where the group velocity becomes extremely superluminal, we can see that δ | |  1 a , which indicates that in these regions the pulse propagation is very much affected by the group dispersion and the group velocity completely looses its physical meaning, even if understood in a purely kinematic sense.