E ﬀ ective Electrodynamics Theory for the Hyperbolic Metamaterial Consisting of Metal–Dielectric Layers

: In this work, we study the dynamical behaviors of the electromagnetic ﬁelds and material responses in the hyperbolic metamaterial consisting of periodically arranged metallic and dielectric layers. The thickness of each unit cell is assumed to be much smaller than the wavelength of the electromagnetic waves, so the e ﬀ ective medium concept can be applied. When electromagnetic (EM) ﬁelds are present, the responses of the medium in the directions parallel to and perpendicular to the layers are similar to those of Drude and Lorentz media, respectively. We derive the time-dependent energy density of the EM ﬁelds and the power loss in the e ﬀ ective medium based on Poynting theorem and the dynamical equations of the polarization ﬁeld. The time-averaged energy density for harmonic ﬁelds was obtained by averaging the energy density in one period, and it reduces to the standard result for the lossless dispersive medium when we turn o ﬀ the loss. A numerical example is given to reveal the general characteristics of the direction-dependent energy storage capacity of the medium. We also show that the Lagrangian density of the system can be constructed. The Euler–Lagrange equations yield the correct dynamical equations of the electromagnetic ﬁelds and the polarization ﬁeld in the medium. The canonical momentum conjugates to every dynamical ﬁeld can be derived from the Lagrangian density via di ﬀ erentiation or variation with respect to that ﬁeld. We apply Legendre transformation to this system and ﬁnd that the resultant Hamiltonian density is identical to the energy density up to an irrelevant divergence term. This coincidence implies the correctness of the energy density formula we obtained before. We also give a brief discussion about the Hamiltonian dynamics description of the system. The Lagrangian description and Hamiltonian formulation presented in this paper can be further developed for studying the elementary excitations or quasiparticles in other hyperbolic metamaterials.


Introduction
Metamaterials usually refer to artificially engineered structures for realizing various unusual optical/electromagnetic properties such as negative refraction [1,2], subwavelength imaging [3], indefinite permittivity [4], near-perfect absorption [5], or invisibility [6,7]. These unusual properties are mainly achieved through the resonance, conductivity, and directionality of the structural components such as split-ring resonators (SRRs), metallic rods array, or subwavelength dielectric-metal multilayers [8]. The resonance and directionality of the constituent components imply that the metamaterials are inherently dispersive, absorptive, and anisotropic. A fundamental problem concerning dispersive media is how to calculate the stored electromagnetic energy density [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. For dispersive media with negligible absorption, the time-averaged energy density as a function of the (complex valued) electric and magnetic fields (in the frequency domain) can be derived by considering the adiabatically varying electromagnetic field [29]. However, such analysis does not work variables. However, up to now, such research in the literature is still lacking. Finally, in the future, more suitable materials may be found and used to design almost lossless hyperbolic metamaterials. If we want to study its quantum behavior, we need to develop a quantum description for the system. However, it is well known that it is only after we formulate the Lagrangian or Hamiltonian of the system correctly that we can quantize the system through some well-established quantization procedure such as Feynman's path integral or the canonical quantization method. The effective field theory developed in this paper may provide useful knowledge for this purpose.

The Effective Fields, Energy Density, and Power Loss
In the first subsection, we will obtain the effective fields in the metal-dielectric multilayer structure according to the boundary conditions of the EM fields and the polarization field (the P field) under the long wavelength approximation. Then, we derive the energy density and the power loss rate in the second subsection using the ED approach based on Poynting's theorem. We will also discuss some important quantities related to monochromatic (single frequency) waves, including the effective permittivity and the time-averaged energy density of the effective medium.
A very important quantity we will encounter in this section is the Q field, which is a newly defined field quantity that appears very naturally through the homogenization process. It replaces the commonly used P field to represent the dynamic part of the D field (the electric displacement field) that does not respond to the E field instantaneously. In other words, in this paper, instead of splitting the D field into D = ε 0 E + P, we split it into D = D ∞ + Q. Here, D ∞ = ε 0 E + P ∞ includes both the ε 0 E term and the instantaneous response P ∞ of the medium. The subscript ∞ means that even at very high frequency, the instantaneous terms still exist. This way of splitting will be found useful when we discuss the derivation of the energy density in the second subsection.

Boundary Conditions and the Effective Fields
We consider a periodic multilayer structure consisting of dielectric and metal layers. Each unit cell has one dielectric and one metal layer (see Figure 1). The lattice constant (the thickness of the unit cell) is a = a m + a d , where a m = f a is the thickness of the metal layer, a d = (1 − f )a is the thickness of the dielectric layer, and f is the filling fraction of the metal layer in one unit cell satisfying 0 < f < 1. The dielectric is a nondispersive material of permittivity ε d = ε 0 ε d , whereas the metal has the Drude-type permittivity ε m (ω) = ε 0 ε m = ε 0 1 − ω 2 p0 ω 2 +iΓω at frequency ω. Here, ε d = ε d /ε 0 and ε m = ε m /ε 0 are the relative (dimensionless) permittivities of the dielectric and metal layers. To make the effective medium theory reasonable, the operating wavelength should be much longer than the lattice constant a. Under this assumption, any field quantity through one single layer does not change value along the direction normal to the layer. The boundary conditions for the E, D, H, and B fields at the dielectric-metal interfaces are the continuity of the tangential component of the E and H (E t ,H t ) fields and the continuity of the normal component of the D and B (D n ,B n ) fields. These boundary conditions are derivable from Maxwell's equations (without external sources) by applying the Stokes and divergence theorem, respectively. We also assume that both the dielectric and metal are non-magnetic materials, so the permeability through the whole structure takes the same value µ 0 as in empty space, and the B field is related to the H field by the simple relation B = µ 0 H. The constitutive relation D = ε 0 E + P between the displacement field D, the electric field E, and the polarization field P is assumed for a single layer as well as for the effective fields in the medium. However, since the boundary conditions for the E and D fields are of different types, the tangential and normal components of the effective fields will be evaluated separately. Hereafter, we denote the direction parallel to and normal to the layers as x and z , respectively. The effective polarization field is the averaged dipole density, which is given by where the superscript "d" and "m" denote the corresponding medium. Similarly, the effective E field is defined by We now define two coefficients x α , z α and a new field Hereafter, we denote the direction parallel to and normal to the layers as x and z, respectively. The effective polarization field is the averaged dipole density, which is given by where the superscript "d" and "m" denote the corresponding medium. Similarly, the effective E field is defined by The boundary conditions for E t and D n lead to the relations In addition, the relation From Equation (1) to Equation (4) and the relation D = ε 0 E + P, we find Similarly, the z component of the D field can be derived as Crystals 2020, 10, 863

of 13
We now define two coefficients α x , α z and a new field and Using these notations, the effective D field in Equations (5) and (6) can be expressed as We have not yet analyzed the dynamical behavior of the P m field. The dynamical behavior of the P m field will determine the dynamical behavior of the effective medium. The Drude-type permittivity p0 ω 2 +iΓω of the metal implies the equation of motion for P m : Taking the x component of Equations (3), (8), and (10), we get the dynamical equation for Q x . .
Here, the effective plasma frequency ω p of the effective medium is defined by the relation Similarly, taking the z component of Equation (10), we get the following dynamical equation However, the right-hand side of Equation (13) must be replaced by the effective fields. This can be done by noting that Equations (2), (3) and (9) tell us Substituting Equation (14) into Equation (13), we get ..
where the resonance frequency ω 0 and the factor F are defined as We are now ready to derive the energy density of the effective medium system. Before doing so, let us check what these equations tell us about the principal permittivity ε x (ω) and ε z (ω). Considering harmonic fields of frequency ω, and replacing all fields with their complex vector representations with time factor e −iωt , under this consideration, Equations (11) and (15) yield Crystals 2020, 10, 863 6 of 13 Here Q x , Q z , E x , and E z are the complex representations of Q x , Q z , E x , and E z . Substituting Equation (17) into Equation (9), and applying the relations D x = ε x E x and D z = ε z E z , we find These results are exactly the same as those obtained directly by using the effective permittivity formula ε at the direction parallel and normal to the layers, as can be easily checked. According to Equations (17) and (18), the Q field tends to zero as the frequency becomes higher and higher. In fact, the Q field is the dynamical part of the P field that does not react immediately to the change of the E field.

Poynting Theorem, Energy Density, and Power Loss
To derive the energy density, we first derive from Maxwell's equations the following equations Now, if the right-hand side of Equation (19) can be written as ∂W ∂t + P loss , and P loss can be identified as the power loss density, then W = W e + W b is the desired energy density of the system, and Equation (19) represents the Poynting theorem (the energy conservation law). Using the simple relation B = µ 0 H, we find H · ∂B ∂t = ∂ ∂t µ 0 2 H 2 , so the magnetic energy density is W b = µ 0 2 H 2 . Furthermore, using Equations (9), (11) and (15), we get Using Equations (8), (12) and (16), we find that the final term of Equation (20) can be re-expressed as a quantity proportional to . P m 2 : Since . P m is the polarization current density, thus Γ is the Joule heat rate density in the metal layer [19,20]. Therefore, it is reasonable to identify the quantity in Equation (21) as the power loss density of the effective medium. With this identification, we can identify the electric energy density W e as Thus, the total energy density is thus given by This result indicates that the total energy density of the system is definitely positive and is consisting of two parts: the non-dispersive part and the dispersive part The different dispersion characters are caused by the facts: W EH is the contribution from the electric-magnetic fields themselves and the part of the polarization field that follows the change of the fields immediately, whereas W Q is originated from the material response that does not follow the fields immediately, because a conduction electron in the metal has nonzero inertial mass.
For harmonic fields, the time average of a product a(t)b(t) (energy density or Poynting vector) can be evaluated by using the formula [37] where T = 2π/ω is the oscillation period, while a and b are the complex representations of the fields a(t) = Re ae −iωt and b(t) = Re be −iωt . Applying Equation (26) to the energy density and power loss density, we get and The second line of Equation (28) is consistent with the generally accepted concept that it is the imaginary part of the permittivity who corresponds to the energy loss in the absorptive medium. Another interesting observation is that when we "turn off" the absorption (i.e., Γ → 0 ), the time averaged energy density becomes The second line of Equation (29) is the prediction to a dispersive medium with negligible absorption, which can be derived by considering the adiabatic variation of the field amplitudes [29].

Numerical Results
We now consider some numerical results of Equations (18), (27) and (28) (see Figure 2). In our simulation, the filling fraction is f = 0.5, and the relative permittivity of the dielectric layers is ε d = 2.5. According to Equations (7), (12) and (16), we find α x = 1.75, α z = 1.4286, F = 2.5, and ω 0 = ω p = 0.5345ω p0 . In addition, we assume that Γ = 0.01ω p0 or equivalently Γ = 0.0187ω p . Substituting these parameters into Equation (17), we get the results in Figure 2a. two frequency ranges are low, although the directionality factor indicates that the absorption property of the medium in these two ranges is also anisotropic. On the other hand, the electric field should be applied to the direction normal to the layers so as to store more energy in the medium if the frequency is within the range between these two specific frequencies. However, for a frequency too close to the resonance frequency 0 ω ω = p , the power loss rate (see Equation (28)) gets so high; thus, the stored energy becomes heat in a short time.  (18)). The green and blue curves are the real parts, whereas the black dotted curve and the broken red curve are their imaginary parts. In (b), we show the directionality factors for the energy density (the blue curve) and power loss (the red brokenline curve). These two factors are defined as the natural log of the ratio between the coefficients of

The Effective Lagrangian Density and Hamiltonian Density
In this section, we will show that within the framework of the Lagrangian-Hamiltonian field theory, the energy density in the effective medium can be identified as the Hamiltonian density of the system. In the first subsection, we will propose a Lagrangian density and show that the EM field  (18) and (27). The relative permittivity of the dielectric layers is assumed to be ε d = 2.5, and the filling fraction is chosen as f = 0.5. In (a), we plot the two principal values of the effective permittivity tensor as functions of frequency. They correspond to the direction parallel and normal to the layers (See Equation (18)). The green and blue curves are the real parts, whereas the black dotted curve and the broken red curve are their imaginary parts. In (b), we show the directionality factors for the energy density (the blue curve) and power loss (the red broken-line curve).
These two factors are defined as the natural log of the ratio between the coefficients of | E z | 2 and | E x | 2 in Equations (27) and (28).
To study the corresponding dispersion/anisotropy and absorption effects in the energy density and power loss, we define two "directionality factors". The first directionality factor is defined as p ω 2 +Γ 2 and β z = α z ε 0 4 1 + and | E z | 2 in Equation (27), respectively. Similarly, the second directionality factor DF loss = ln γ z γ x is defined according to Equation (28), where γ x = α x Γε 0 2 ω 2 p ω 2 +Γ 2 = ω 2 Im(ε x ) and These two factors as functions of frequency are shown in Figure 2b. According to these numerical results, the medium at an angular frequency lower than 0.8ω p or higher than 1.2ω p can be treated as a hyperbolic metamaterial because the real parts of the ε x and ε z have different signs, and the imaginary parts of them are small enough to be negligible. In the frequency range close to ω p , the imaginary parts become so large that it is not suitable to treat this effective medium as a hyperbolic metamaterial, although our energy density formula can still be applied if the long wavelength approximation is accurate enough. Furthermore, according to the simulations of the directionality factors shown in Figure 2b, the energy storage ability of the effective medium becomes isotropic with respect to the electric field at two specific angular frequencies: 0.489ω p and 3.124ω p . Below 0.489ω p or above 3.124ω p , the medium can store more energy if the electric field is applied in the directions parallel to the layers. The power loss corresponding to these two frequency ranges are low, although the directionality factor indicates that the absorption property of the medium in these two ranges is also anisotropic. On the other hand, the electric field should be applied to the direction normal to the layers so as to store more energy in the medium if the frequency is within the range between these two specific frequencies. However, for a frequency too close to the resonance frequency ω p = ω 0 , the power loss rate (see Equation (28)) gets so high; thus, the stored energy becomes heat in a short time.

The Effective Lagrangian Density and Hamiltonian Density
In this section, we will show that within the framework of the Lagrangian-Hamiltonian field theory, the energy density in the effective medium can be identified as the Hamiltonian density of the system. In the first subsection, we will propose a Lagrangian density and show that the EM field equations and the dynamic equations of the Q field can be obtained from the Euler-Lagrange equations. Since the medium is absorptive, here, we use the modified version of the Euler-Lagrange equations, which include not only the Lagrangian density itself but also Rayleigh's dissipation function density [36]. Then, the Hamiltonian density is obtained in the second subsection by applying the standard Legendre transformation. Here, we must emphasize that in this framework, the frequently used E, H, D, and B fields are not the dynamic variables but derived fields. The proper dynamic variables are ψ α = [ϕ, A x , A z , Q x , Q z ]; here, ϕ and A are the scalar and vector potential, respectively.

Lagrangian Density and Euler-Lagrange Equations
In this section, we will construct the Lagrangian density for the effective medium system as a function of the scalar potential ϕ, the vector potential A, the Q field, and their time and space derivatives. It is clear that in Maxwell's equations, the Gauss law ∇ · B = 0 for the B field and the Faraday's induction law ∇ × E = −∂B/∂t are automatically satisfied because B = ∇ × A implies ∇ · B = 0 and E = −∇ϕ − . A implies ∇ × E = −∂B/∂t. We will show that the other two Maxwell's equations ∇ · D = 0 and ∇ × H = ∂D/∂t as well as the equations of motion for the Q fields (Equations (11) and (15)) can be derived from the Euler-Lagrange equations (with dissipation) [36] Here, ψ α = [ϕ, A x , A z , Q x , Q z ] are the dynamical fields involved in the Lagrangian density, and the relations D x = α x ε 0 E x + Q x and D z = α z ε 0 E z + Q z were used. The dissipation function density F in Equation (30) is defined as which has the value equal to one-half of the power loss. Observing Equations (24) and (25), it is not difficult to guess the Lagrangian density. It should be a sum of three kinds of terms where L EH describes the electromagnetic fields and L C describes the matter-field coupling As mentioned before, using this Lagrangian density L and the dissipation function density F in Equation (31), one can derive all the dynamical equations of the effective fields from Equation (30). It can be checked that for ψ α = ϕ we get the Gauss law ∇ · D = 0; for ψ α = A , we get the Ampere's law ∇ × H = ∂D/∂t; and for ψ α = Q , we get the equations of motion for the Q field (i.e., Equations (11) and (15)).

Canonical Momenta, Legendre Transformation, and Hamiltonian Density
In order to derive the Hamiltonian density, we have to derive all the canonical momenta first. The canonical momentum π α conjugate to the dynamical field ψ α is defined by π α = ∂L ∂ . ψ α . For of the dielectric layers is assumed to be larger than 1 even at high frequency. This assumption also leads to the necessity of splitting the effective displacement field D as D = D ∞ + Q instead of D = ε 0 E + P, and the prediction about the directional dependent behavior in the energy storage ability that was discussed in the previous section.

Conclusions
In this paper, we have completed the derivations of the effective fields, energy density, Lagrangian density, and Hamiltonian density for the electrodynamics in the effective medium that consists of dielectric-metal layers. We have also discussed how to obtain the frequency domain quantities such the permittivities and time-averaged energy density from our time domain formulas. It is found that the Hamiltonian density is the same as the energy density, up to an irrelevant divergence term. This fact confirms the correctness of the energy density formula. One important finding in this study is that the dynamical part of the electric displacement field-the Q field-is a more convenient dynamical field than the usual P field (the polarization field) when studying the electrodynamics of this effective medium.
The Lagrangian field theory is a systematic method that yields the correct equations of motion for the polarization field and the vector potential through the Euler-Lagrange equations. Since the system is dissipative, a dissipation function was introduced that takes care of the effect related to energy loss. The Lagrangian/Hamiltonian field theory framework developed here may provide the essential knowledge for further developing the quantum theory of this kind of media.