Effective Electrodynamics Theory for the Hyperbolic Metamaterial Consisting of Metal–Dielectric Layers

Abstract

In this work, we study the dynamical behaviors of the electromagnetic fields 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 effective medium concept can be applied. When electromagnetic (EM) fields 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 fields and the power loss in the effective medium based on Poynting theorem and the dynamical equations of the polarization field. The time-averaged energy density for harmonic fields 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 off 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 fields and the polarization field in the medium. The canonical momentum conjugates to every dynamical field can be derived from the Lagrangian density via differentiation or variation with respect to that field. We apply Legendre transformation to this system and find 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.

1. 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 when finite absorption is present. Recently, two different approaches were proposed to resolve such a non-trivial problem [16,18]. For the wire–SRR metamaterials, time-averaged energy density formula can be derived using the equivalent circuit (EC) method [10,16,17]. On the other hand, instantaneous energy density formula (the time domain formula) can be obtained using the electrodynamics (ED) method with Poynting theorem [9,12,13,14,15,18,19,20,24,27,28]. However, quite a few controversies exist in the literature, and they need to be resolved [16,17,18,19,20]. In our previous studies, we found that the time domain formula for energy density can be uniquely determined by the ED approach provided we know how to identify the power loss [19,20]. Recently, we also developed the Lagrangian field theory description for the wire–SRR and chiral metamaterials [30]. In this framework, the Hamiltonian density for the metamaterial can be obtained through the Legendre transformation. It is found that the Hamiltonian densities for the dynamical fields in these metamaterials are the same as the energy densities we already obtained, up to some irrelevant divergence terms.

Recently, hyperbolic metamaterials belonging to the category of anisotropic media have received more and more attention [31,32,33,34,35]. The hyperbolic dispersion of this kind of metamaterial is caused by the opposite signs of the two principal values of the permittivity tensor along and perpendicular to the optic axis. One of the most important applications of hyperbolic metamaterial is the hyperlens [34], which can image subwavelength objects in the far-field region, overcoming the shortcomings of the superlens [1,3] that can only image the same objects in the near-field region. The simplest example of a hyperbolic metamaterial is a dielectric–metal multilayer structure that operates under the long-wavelength limit [35]. Although it seems that such metamaterials are much simpler than the wire–SRR and chiral metamaterials, to the best of our knowledge, the energy density problem and the Lagrangian description for this kind of metamaterial have not yet been studied. This is the main motivation of our present study. In addition, a better understanding about the energy density, Lagrangian description, and Hamiltonain theory for the hyperbolic metamaterial can help researchers further explore the dynamical behaviors of EM waves in this kind of media without being restricted to the frequency domain phenomena. This may have practical importance in the future. The Lagrangian and Hamiltonian descriptions also provide a starting point for the development of quantum description of the electrodynamics in the metamaterials [30].

In this paper, we study the effective electrodynamics of dielectric–metal multilayer structures under the limitation of long wavelengths. Assuming that each metal layer is absorptive, we study the energy density and power loss of the system. We first discuss the boundary conditions of the dynamic fields in each dielectric layer and metal layer, and we derive the effective fields. We also derive the effective permittivities in the frequency domain based on our theory of the effective fields. Then, we study the dynamical evolution equations of the effective electric, magnetic, displacement, and polarization fields and derive the effective energy density using the ED method. We discuss the loss effect and provide a suitable dissipation function [36] so that we can derive the correct equations of motion for the dynamic fields from the Euler–Lagrange equation with dissipation. The Hamiltonian density is derived by applying the Legendre transformation to this system. The resultant Hamiltonian density is found to be the same as the energy density obtained before, up to a divergence term.

The Lagrangian description for the electrodynamics of metamaterials has many advantages. First of all, the Lagrangian description provides a well-established approach—the Legendre transformation, for us to derive the Hamiltonian density, which is equal to the desired energy density we are looking for. Second, the conservation law related to the symmetry of the Lagrangian can be easily obtained by applying Noether’s theorem [36]. For example, the Poynting theorem of lossless electrodynamic systems is a consequence of the translation symmetry of the system in space time. In addition, Noether’s theorem can also provide us with the correct formula for the stress tensor of the fields (we will not get the formula in this paper) [29], which can help us analyze the momentum transfer caused by light–matter interaction. Third, the use of the Lagrangian description in conjunction with Rayleigh’s dissipation function also provides a convenient way to study the dynamical phenomena in absorptive systems whose power loss is proportional to the square of some velocity or current of the dynamical 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.

2. 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={\epsilon }_{0}E+P$, we split it into $D=D_{\infty }+Q$. Here, $D_{\infty }={\epsilon }_{0}E+P_{\infty }$ includes both the ${\epsilon }_{0}E$ term and the instantaneous response $P_{\infty }$ of the medium. The subscript $\infty$ 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.

2.1. 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=fa$ 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 ${\epsilon }_d={\epsilon }_0{\tilde{\epsilon }}_d$, whereas the metal has the Drude-type permittivity ${\epsilon }_m(\omega )={\epsilon }_0{\tilde{\epsilon }}_m={\epsilon }_0\left(1-\frac{{\omega }_{p0}^2}{{\omega }^2+i\mathsf{\Gamma }\omega }\right)$ at frequency $\omega$. Here, ${\tilde{\epsilon }}_d={\epsilon }_d/{\epsilon }_0$ and ${\tilde{\epsilon }}_m={\epsilon }_m/{\epsilon }_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 ${\mu }_0$ as in empty space, and the $B$ field is related to the $H$ field by the simple relation $B={\mu }_{0}H$. The constitutive relation $D={\epsilon }_{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

$$P=\left(1-f\right)P^d+f{P}^m,$$

(1)

where the superscript “d” and “m” denote the corresponding medium. Similarly, the effective $E$ field is defined by

$$E=\left(1-f\right)E^d+f{E}^m.$$

(2)

The boundary conditions for $E_t$ and $D_n$ lead to the relations

$$E_x=E_x^d=E_x^m=D_x^d/{\epsilon }_d,D_z=D_z^d=D_z^m={\epsilon }_d{E}_z^d.$$

(3)

In addition, the relation $D^d={\epsilon }_0{E}^d+P^d={\epsilon }_d{E}^d$ implies

$$P_x^d=\left({\tilde{\epsilon }}_d-1\right){\epsilon }_0{E}_x,P_z^d=\left(1-\frac{1}{{\tilde{\epsilon }}_d}\right)D_z.$$

(4)

From Equation (1) to Equation (4) and the relation $D={\epsilon }_{0}E+P$, we find

$$\begin{array}{ll}{D}_x& ={\epsilon }_0{E}_x+\left(1-f\right)P_x^d+f{P}_x^m\\ & ={\epsilon }_0{E}_x+\left(1-f\right)\left({\tilde{\epsilon }}_d-1\right){\epsilon }_0{E}_x+f{P}_x^m.\\ & ={\epsilon }_0\left[f+\left(1-f\right){\tilde{\epsilon }}_d\right]E_x+f{P}_x^m\end{array}$$

(5)

Similarly, the z component of the $D$ field can be derived as

$$\begin{array}{ll}{D}_z& ={\epsilon }_0{E}_z+\left(1-f\right)P_z^d+f{P}_z^m\\ & ={\epsilon }_0{E}_z+\left(1-f\right)\left(D_z-{\epsilon }_0{E}_z\right)+f{P}_x^m\\ & ={\epsilon }_0{E}_z+\left(1-f\right)\left(1-\frac{1}{{\tilde{\epsilon }}_d}\right)D_z+f{P}_z^m\\ & =\frac{{\tilde{\epsilon }}_d\left({\epsilon }_0{E}_z+f{P}_z^m\right)}{f{\tilde{\epsilon }}_d+\left(1-f\right)}\end{array}.$$

(6)

We now define two coefficients ${\alpha }_x$, ${\alpha }_z$ and a new field $Q=[Q_x,Q_z]$ as

$${\alpha }_x=f+(1-f){\tilde{\epsilon }}_d,{\alpha }_z=\frac{{\tilde{\epsilon }}_d}{f{\tilde{\epsilon }}_d+(1-f)}$$

(7)

and

$$Q_x=f{P}_x^m,Q_z={\alpha }_{z}f{P}_z^m.$$

(8)

Using these notations, the effective $D$ field in Equations (5) and (6) can be expressed as

$$D_x={\alpha }_x{\epsilon }_0{E}_x+Q_x,D_z={\alpha }_z\left({\epsilon }_0{E}_z+f{P}_z^m\right)={\alpha }_z{\epsilon }_0{E}_z+Q_z.$$

(9)

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 ${\epsilon }_m(\omega )={\epsilon }_0\left(1-\frac{{\omega }_{p0}^2}{{\omega }^2+i\mathsf{\Gamma }\omega }\right)$ of the metal implies the equation of motion for $P^m$:

$${\ddot{\mathbf{P}}}^m+{\Gamma }{\dot{\mathbf{P}}}^m={\epsilon }_0{\omega }_{p0}^2{E}^m.$$

(10)

Taking the x component of Equations (3), (8), and (10), we get the dynamical equation for $Q_x$

$${\ddot{Q}}_x+{\Gamma }{\dot{Q}}_x={\epsilon }_{0}f{\omega }_{p0}^2{E}_x={\alpha }_x{\epsilon }_0{\omega }_p^2{E}_x.$$

(11)

Here, the effective plasma frequency ${\omega }_p$ of the effective medium is defined by the relation

$${\omega }_p^2=\frac{f{\omega }_{p0}^2}{{\alpha }_x}.$$

(12)

Similarly, taking the z component of Equation (10), we get the following dynamical equation

$${\ddot{Q}}_z+{\Gamma }{\dot{Q}}_z={\alpha }_z{\epsilon }_0{\omega }_{p0}^{2}f{E}_z^m.$$

(13)

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

$$f{E}_z^m=E_z-(1-f)\frac{D_z}{{\epsilon }_d}={\alpha }_{z}f{E}_z-\frac{(1-f)}{{\epsilon }_d}{Q}_z.$$

(14)

Substituting Equation (14) into Equation (13), we get

$${\ddot{Q}}_z+{\Gamma }{\dot{Q}}_z+{\omega }_0^2{Q}_z=F{\alpha }_z{\epsilon }_0{\omega }_0^2{E}_z,$$

(15)

where the resonance frequency ${\omega }_0$ and the factor $F$ are defined as

$${\omega }_0^2=\frac{(1-f){\alpha }_z}{{\tilde{\epsilon }}_d}{\omega }_{p0}^2,F=\frac{{\omega }_{p0}^2}{{\omega }_0^2}-1=\frac{f{\tilde{\epsilon }}_d}{1-f}.$$

(16)

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 ${\epsilon }_x(\omega )$ and ${\epsilon }_z(\omega )$. Considering harmonic fields of frequency $\omega$, and replacing all fields with their complex vector representations with time factor $e^{-i\omega t}$, under this consideration, Equations (11) and (15) yield

$${\tilde{Q}}_x=-\frac{{\alpha }_x{\epsilon }_0{\omega }_p^2}{{\omega }^2+i{\Gamma }\omega }{\tilde{E}}_x,{\tilde{Q}}_z=-\frac{F{\alpha }_z{\epsilon }_0{\omega }_0^2}{{\omega }^2-{\omega }_0^2+i{\Gamma }\omega }{\tilde{E}}_z(\mathrm{complex}\mathrm{rep}.).$$

(17)

Here ${\tilde{Q}}_x$, ${\tilde{Q}}_z$, ${\tilde{E}}_x$, and ${\tilde{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 ${\tilde{D}}_x={\epsilon }_x{\tilde{E}}_x$ and ${\tilde{D}}_z={\epsilon }_z{\tilde{E}}_z$, we find

$${\epsilon }_x={\alpha }_x{\epsilon }_0\left(1-\frac{{\omega }_p^2}{{\omega }^2+i{\Gamma }\omega }\right),{\epsilon }_z={\alpha }_z{\epsilon }_0\left(1-\frac{F{\omega }_0^2}{{\omega }^2-{\omega }_0^2+i{\Gamma }\omega }\right).$$

(18)

These results are exactly the same as those obtained directly by using the effective permittivity formula ${\epsilon }_x=f{\epsilon }_m+(1-f){\epsilon }_d$ and $1/{\epsilon }_z=f/{\epsilon }_m+(1-f)/{\epsilon }_d$ 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.

2.2. Poynting Theorem, Energy Density, and Power Loss

To derive the energy density, we first derive from Maxwell’s equations the following equations

$$-\nabla \cdot S=E\cdot \frac{\partial D}{\partial t}+H\cdot \frac{\partial B}{\partial t}.$$

(19)

Now, if the right-hand side of Equation (19) can be written as $\frac{\partial W}{\partial 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={\mu }_{0}H$, we find $H\cdot \frac{\partial B}{\partial t}=\frac{\partial }{\partial t}\left(\frac{{\mu }_0}{2}{H}^2\right)$, so the magnetic energy density is $W_b=\frac{{\mu }_0}{2}{H}^2$. Furthermore, using Equations (9), (11) and (15), we get

$$\begin{array}{ll}E\cdot \frac{\partial D}{\partial t}& =E_x\frac{\partial D_x}{\partial t}+E_z\frac{\partial D_z}{\partial t}=E_x\frac{\partial }{\partial t}\left({\alpha }_x{\epsilon }_0{E}_x+Q_x\right)+E_z\frac{\partial }{\partial t}\left({\alpha }_z{\epsilon }_0{E}_z+Q_z\right)\\ & =\frac{\partial }{\partial t}\left[\frac{{\epsilon }_0}{2}\left({\alpha }_x{E}_x^2+{\alpha }_z{E}_z^2\right)\right]+\frac{\left({\ddot{Q}}_x+\Gamma {\dot{Q}}_x\right){\dot{Q}}_x}{{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{\left({\ddot{Q}}_z+\Gamma {\dot{Q}}_z+{\omega }_0^2{Q}_z\right){\dot{Q}}_z}{F{\alpha }_z{\epsilon }_0{\omega }_0^2}\\ & =\frac{\partial }{\partial t}\left[\frac{{\epsilon }_0}{2}\left({\alpha }_x{E}_x^2+{\alpha }_z{E}_z^2\right)+\frac{{\dot{Q}}_x^2}{2{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{{\dot{Q}}_z^2+{\omega }_0^2{Q}_z^2}{2F{\alpha }_z{\epsilon }_0{\omega }_0^2}\right]+{\Gamma }\left(\frac{{\dot{Q}}_x^2}{{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{{\dot{Q}}_z^2}{F{\alpha }_z{\epsilon }_0{\omega }_0^2}\right)\end{array}.$$

(20)

Using Equations (8), (12) and (16), we find that the final term of Equation (20) can be re-expressed as a quantity proportional to ${\left({\dot{\mathbf{P}}}^m\right)}^2$:

$$P_{loss}={\Gamma }\left(\frac{{\dot{Q}}_x^2}{{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{{\dot{Q}}_z^2}{F{\alpha }_z{\epsilon }_0{\omega }_0^2}\right)=\frac{f{\Gamma }}{{\epsilon }_0{\omega }_{p0}^2}{\left({\dot{\mathbf{P}}}^m\right)}^2.$$

(21)

Since ${\dot{\mathbf{P}}}^m$ is the polarization current density, thus $\frac{{\Gamma }}{{\epsilon }_0{\omega }_{p0}^2}{\left({\dot{\mathbf{P}}}^m\right)}^2$ 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

$$W_e=\frac{{\epsilon }_0}{2}\left({\alpha }_x{E}_x^2+{\alpha }_z{E}_z^2\right)+\frac{{\dot{Q}}_x^2}{2{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{{\dot{Q}}_z^2+{\omega }_0^2{Q}_z^2}{2F{\alpha }_z{\epsilon }_0{\omega }_0^2}.$$

(22)

Thus, the total energy density is thus given by

$$W=W_e+W_b=\frac{{\epsilon }_0}{2}\left({\alpha }_x{E}_x^2+{\alpha }_z{E}_z^2\right)+\frac{{\mu }_0}{2}{H}^2+\frac{{\dot{Q}}_x^2}{2{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{{\dot{Q}}_z^2+{\omega }_0^2{Q}_z^2}{2F{\alpha }_z{\epsilon }_0{\omega }_0^2}.$$

(23)

This result indicates that the total energy density of the system is definitely positive and is consisting of two parts: the non-dispersive part

$$W_{EH}=\frac{{\epsilon }_0}{2}\left({\alpha }_x{E}_x^2+{\alpha }_z{E}_z^2\right)+\frac{{\mu }_0}{2}{H}^2,$$

(24)

and the dispersive part

$$W_Q=\frac{{\dot{Q}}_x^2}{2{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{{\dot{Q}}_z^2+{\omega }_0^2{Q}_z^2}{2F{\alpha }_z{\epsilon }_0{\omega }_0^2}.$$

(25)

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]

$$\langle a(t)b(t)\rangle =\frac{1}{T}{\displaystyle {\int }_0^{T}a(t)b(t)dt=}\frac{1}{2}\mathrm{Re}\left(\tilde{a}{\tilde{b}}^{*}\right),$$

(26)

where $T=2\pi /\omega$ is the oscillation period, while $\tilde{a}$ and $\tilde{b}$ are the complex representations of the fields $a(t)=\mathrm{Re}\left(\tilde{a}{e}^{-i\omega t}\right)$ and $b(t)=\mathrm{Re}\left(\tilde{b}{e}^{-i\omega t}\right)$.

Applying Equation (26) to the energy density and power loss density, we get

$$\begin{array}{ll}\langle W\rangle & =\frac{{\epsilon }_0}{4}\left({\alpha }_x|{\tilde{E}}_x{|}^2+{\alpha }_z|{\tilde{E}}_z{|}^2\right)+\frac{{\mu }_0}{4}|\tilde{H}{|}^2+\frac{{\omega }^2|{\tilde{Q}}_x{|}^2}{4{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{\left({\omega }^2+{\omega }_0^2\right)|{\tilde{Q}}_z{|}^2}{4F{\alpha }_z{\epsilon }_0{\omega }_0^2}\\ & =\frac{{\alpha }_x{\epsilon }_0}{4}\left(1+\frac{{\omega }_p^2}{{\omega }^2+{\mathsf{\Gamma }}^2}\right)|{\tilde{E}}_x{|}^2+\frac{{\alpha }_z{\epsilon }_0}{4}\left[1+\frac{F{\omega }_0^2\left({\omega }^2+{\omega }_0^2\right)}{{\left({\omega }^2-{\omega }_0^2\right)}^2+{\mathsf{\Gamma }}^2{\omega }^2}\right]|{\tilde{E}}_z{|}^2+\frac{{\mu }_0}{4}|\tilde{H}{|}^2\end{array}$$

(27)

and

$$\begin{array}{ll}\langle P_{loss}\rangle & =\frac{\mathsf{\Gamma }{\epsilon }_0}{2}\left[\frac{{\alpha }_x{\omega }_p^2}{{\omega }^2+{\mathsf{\Gamma }}^2}|{\tilde{E}}_x{|}^2+\frac{{\alpha }_{z}F{\omega }_0^2{\omega }^2}{{\left({\omega }^2-{\omega }_0^2\right)}^2+{\mathsf{\Gamma }}^2{\omega }^2}|{\tilde{E}}_z{|}^2\right]\\ & =\frac{\omega }{2}\mathrm{Im}\left({\epsilon }_x|{\tilde{E}}_x{|}^2+{\epsilon }_z|{\tilde{E}}_z{|}^2\right)\end{array}.$$

(28)

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., $\mathsf{\Gamma }\to 0$), the time averaged energy density becomes

$$\begin{array}{ll}{\langle W\rangle }_{\mathsf{\Gamma }\to 0}& =\frac{{\alpha }_x{\epsilon }_0}{4}\left(1+\frac{{\omega }_p^2}{{\omega }^2}\right)|{\tilde{E}}_x{|}^2+\frac{{\alpha }_z{\epsilon }_0}{4}\left[1+\frac{F{\omega }_0^2\left({\omega }^2+{\omega }_0^2\right)}{{\left({\omega }^2-{\omega }_0^2\right)}^2}\right]|{\tilde{E}}_z{|}^2+\frac{{\mu }_0}{4}|\tilde{H}{|}^2\\ & =\frac{1}{4}\left[\frac{\partial {\left(\omega {\epsilon }_x\right)}_{\mathsf{\Gamma }\to 0}}{\partial \omega }|{\tilde{E}}_x{|}^2+\frac{\partial {\left(\omega {\epsilon }_z\right)}_{\mathsf{\Gamma }\to 0}}{\partial \omega }|{\tilde{E}}_z{|}^2\right]\end{array}$$

(29)

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].

2.3. 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 ${\tilde{\epsilon }}_d=2.5$. According to Equations (7), (12) and (16), we find ${\alpha }_x=1.75$, ${\alpha }_z=1.4286$, $F=2.5$, and ${\omega }_0={\omega }_p=0.5345{\omega }_{p0}$. In addition, we assume that $\mathsf{\Gamma }=0.01{\omega }_{p0}$ or equivalently $\mathsf{\Gamma }=0.0187{\omega }_p$. Substituting these parameters into Equation (17), we get the results in Figure 2a.

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 $D{F}_{energy}=\ln \left(\frac{{\beta }_z}{{\beta }_x}\right)$, where ${\beta }_x=\frac{{\alpha }_x{\epsilon }_0}{4}\left(1+\frac{{\omega }_p^2}{{\omega }^2+{\mathsf{\Gamma }}^2}\right)$ and ${\beta }_z=\frac{{\alpha }_z{\epsilon }_0}{4}\left[1+\frac{F{\omega }_0^2\left({\omega }^2+{\omega }_0^2\right)}{{\left({\omega }^2-{\omega }_0^2\right)}^2+{\mathsf{\Gamma }}^2{\omega }^2}\right]$ are the coefficients of ${|{\tilde{E}}_x|}^2$ and ${|{\tilde{E}}_z|}^2$ in Equation (27), respectively. Similarly, the second directionality factor $D{F}_{loss}=\ln \left(\frac{{\gamma }_z}{{\gamma }_x}\right)$ is defined according to Equation (28), where ${\gamma }_x=\frac{{\alpha }_x\mathsf{\Gamma }{\epsilon }_0}{2}\left(\frac{{\omega }_p^2}{{\omega }^2+{\mathsf{\Gamma }}^2}\right)=\frac{\omega }{2}\mathrm{Im}\left({\epsilon }_x\right)$ and ${\gamma }_z=\frac{{\alpha }_z\mathsf{\Gamma }{\epsilon }_0}{2}\left[\frac{F{\omega }_0^2{\omega }^2}{{\left({\omega }^2-{\omega }_0^2\right)}^2+{\mathsf{\Gamma }}^2{\omega }^2}\right]=\frac{\omega }{2}\mathrm{Im}\left({\epsilon }_z\right)$ are the coefficients of ${|{\tilde{E}}_x|}^2$ of ${|{\tilde{E}}_z|}^2$. 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{\omega }_p$ or higher than $1.2{\omega }_p$ can be treated as a hyperbolic metamaterial because the real parts of the ${\epsilon }_x$ and ${\epsilon }_z$ have different signs, and the imaginary parts of them are small enough to be negligible. In the frequency range close to ${\omega }_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{\omega }_p$ and $3.124{\omega }_p$. Below $0.489{\omega }_p$ or above $3.124{\omega }_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 ${\omega }_p={\omega }_0$, the power loss rate (see Equation (28)) gets so high; thus, the stored energy becomes heat in a short time.

3. 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 ${\psi }_{\alpha }=[\phi ,A_x,A_z,Q_x,Q_z]$; here, $\phi$ and $A$ are the scalar and vector potential, respectively.

3.1. 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 $\phi$, the vector potential $A$, the $Q$ field, and their time and space derivatives. It is clear that in Maxwell’s equations, the Gauss law $\nabla \cdot B=0$ for the $B$ field and the Faraday’s induction law $\nabla \times E=-\partial B/\partial t$ are automatically satisfied because $B=\nabla \times A$ implies $\nabla \cdot B=0$ and $E=-\nabla \phi -\dot{A}$ implies $\nabla \times E=-\partial B/\partial t$. We will show that the other two Maxwell’s equations $\nabla \cdot D=0$ and $\nabla \times H=\partial D/\partial 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]

$${\partial }_t\left(\frac{\partial \mathcal{L}}{\partial \left({\partial }_t{\psi }_{\alpha }\right)}\right)+{\partial }_j\left(\frac{\partial \mathcal{L}}{\partial \left({\partial }_j{\psi }_{\alpha }\right)}\right)-\frac{\partial \mathcal{L}}{\partial {\psi }_{\alpha }}=-\frac{\partial \mathcal{F}}{\partial \left({\partial }_t{\psi }_{\alpha }\right)}.$$

(30)

Here, ${\psi }_{\alpha }=[\phi ,A_x,A_z,Q_x,Q_z]$ are the dynamical fields involved in the Lagrangian density, and the relations $D_x={\alpha }_x{\epsilon }_0{E}_x+Q_x$ and $D_z={\alpha }_z{\epsilon }_0{E}_z+Q_z$ were used. The dissipation function density $\mathcal{F}$ in Equation (30) is defined as

$$\mathcal{F}=\frac{\mathsf{\Gamma }}{2}\left(\frac{{\dot{Q}}_x^2}{{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{{\dot{Q}}_z^2}{F{\alpha }_z{\epsilon }_0{\omega }_0^2}\right)=\frac{1}{2}{P}_{loss}$$

(31)

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

$$\mathcal{L}={\mathcal{L}}_{EH}+{\mathcal{L}}_Q+{\mathcal{L}}_C,$$

(32)

where ${\mathcal{L}}_{EH}$ describes the electromagnetic fields

$${\mathcal{L}}_{EH}=\frac{{\epsilon }_0}{2}\left[{\alpha }_x{\left({\partial }_x\phi +{\dot{A}}_x\right)}^2+{\alpha }_z{\left({\partial }_z\phi +{\dot{A}}_z\right)}^2\right]-\frac{1}{2{\mu }_0}{\left(\nabla \times A\right)}^2,$$

(33a)

${\mathcal{L}}_Q$ describes the $Q$ field

$${\mathcal{L}}_Q=\frac{{\dot{Q}}_x^2}{2{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{{\dot{Q}}_z^2-{\omega }_0^2{Q}_z^2}{2F{\alpha }_z{\epsilon }_0{\omega }_0^2},$$

(33b)

and ${\mathcal{L}}_C$ describes the matter–field coupling

$${\mathcal{L}}_C=E_x{Q}_x+E_z{Q}_z=-\left({\partial }_x\phi +{\dot{A}}_x\right)Q_x-\left({\partial }_z\phi +{\dot{A}}_z\right)Q_z.$$

(33c)

As mentioned before, using this Lagrangian density $\mathcal{L}$ and the dissipation function density $\mathcal{F}$ in Equation (31), one can derive all the dynamical equations of the effective fields from Equation (30). It can be checked that for ${\psi }_{\alpha }=\phi$ we get the Gauss law $\nabla \cdot D=0$; for ${\psi }_{\alpha }=A$, we get the Ampere’s law $\nabla \times H=\partial D/\partial t$; and for ${\psi }_{\alpha }=Q$, we get the equations of motion for the $Q$ field (i.e., Equations (11) and (15)).

3.2. 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 ${\pi }_{\alpha }$ conjugate to the dynamical field ${\psi }_{\alpha }$ is defined by ${\pi }_{\alpha }=\frac{\partial \mathcal{L}}{\partial {\dot{\psi }}_{\alpha }}$. For ${\psi }_{\alpha }=[\phi ,A_x,A_z,Q_x,Q_z]$, we find ${\pi }_{\alpha }=[0,-D_x,-D_z,{\dot{Q}}_x/({\alpha }_x{\epsilon }_0{\omega }_p^2),{\dot{Q}}_z/({\alpha }_{z}F{\epsilon }_0{\omega }_0^2)]$. It is interesting to note that the canonical momentum ${\pi }_{\phi }$ conjugate to $\phi$ is equal to zero. This is a common feature of electrodynamics systems, which reflects the fact that $\phi$ is a redundant dynamical variable.

The Hamiltonian density $\mathcal{H}$ is given by the Legendre transformation

$$\mathcal{H}={\pi }_{\alpha }{\dot{\psi }}_{\alpha }-\mathcal{L}=\frac{{\epsilon }_0}{2}\left({\alpha }_x{E}_x^2+{\alpha }_z{E}_z^2\right)+\frac{{\mu }_0}{2}{H}^2+\frac{{\dot{Q}}_x^2}{2{\alpha }_x{\epsilon }_0{\omega }_p^2}+\frac{{\dot{Q}}_z^2+{\omega }_0^2{Q}_z^2}{2F{\alpha }_z{\epsilon }_0{\omega }_0^2}+D\cdot \nabla \phi .$$

(34)

According to this result, the Hamiltonian density is almost the same as the energy density we obtained in Equation (23). The only difference is the new term $D\cdot \nabla \phi =\nabla \cdot \left(\phi D\right)-\phi \nabla \cdot D=\nabla \cdot \left(\phi D\right)$. Here, we have used the Gauss law $\nabla \cdot D=0$. This additional divergence term can be dropped from Equation (34), because it is only a “surface term” and will not influence the dynamical equations for the effective fields. After dropping this surface term, the final term in Equation (34) can be replaced by a term of zero value: $-\phi \nabla \cdot D$. In the framework of Dirac’s treatment to the constraint systems [38], we can treat $\phi$ as a Lagrangian multiplier, and $\nabla \cdot D=0$ is the constraint to the canonical momentum ${\pi }_A=-D$.

3.3. Hamiltonian Equations of Motion for the Canonical Fields and Momenta

We now study the problem of how to express the dynamical behaviors of this effective medium system in the framework of Hamiltonian dynamics. This is an important step if we want to quantize this system using the canonical quantization method. Since the canonical momenta conjugate to ${\psi }_{\alpha }=[\phi ,A_x,A_z,Q_x,Q_z]$ are ${\pi }_{\alpha }=[0,-D_x,-D_z,{\dot{Q}}_x/({\alpha }_x{\epsilon }_0{\omega }_p^2),{\dot{Q}}_z/({\alpha }_{z}F{\epsilon }_0{\omega }_0^2)]$, respectively, we can rewrite the Hamiltonian density formula (Equation (34)) as

$$\mathcal{H}=\frac{{\left({\pi }_{A_x}+Q_x\right)}^2}{2{\epsilon }_0{\alpha }_x}+\frac{{\left({\pi }_{A_z}+Q_z\right)}^2}{2{\epsilon }_0{\alpha }_z}+\frac{{\left(\nabla \times A\right)}^2}{2{\mu }_0}+\frac{{\epsilon }_0{\alpha }_x{\omega }_p^2{\pi }_{Q_x}^2}{2}+\frac{{\epsilon }_{0}F{\alpha }_z{\omega }_0^2{\pi }_{Q_z}^2}{2}+\frac{Q_z^2}{2F{\alpha }_z{\epsilon }_0}-{\pi }_A\cdot \nabla \phi .$$

(35)

Now, the total Hamiltonian is defined as the volume integration of the Hamiltonian density

$$H={\displaystyle \int \mathcal{H}}{d}^{3}x,$$

(36)

and the dynamical equations for the ${\psi }_{\alpha }$ and ${\pi }_{\alpha }$ are

$${\dot{\psi }}_{\alpha }=\frac{\delta H}{\delta {\pi }_{\alpha }},{\dot{\pi }}_{\alpha }=-\frac{\delta H}{\delta {\psi }_{\alpha }}-\frac{\partial \mathcal{F}}{\partial {\dot{\psi }}_{\alpha }}.$$

(37)

Here, the notation such as $\frac{\delta H}{\delta {\pi }_{\alpha }}$ denotes the functional derivative [39], whereas $\frac{\partial \mathcal{F}}{\partial {\dot{\psi }}_{\alpha }}$ is simply the usual partial derivative. The appearance of the term $-\frac{\partial \mathcal{F}}{\partial {\dot{\psi }}_{\alpha }}$ is due to the dissipation effect of the system. One can check that Equation (37) indeed reproduces all the dynamical equations for the dynamical fields. If the loss of the system is so small to be negligible, this additional term disappears, and Equation (37) reduces to the standard Hamilton equations, and we can quantize it straightforwardly. However, if the loss effect cannot be ignored, this term must be kept, and we do not really know how to deal with it in the Hamiltonian dynamics framework. We left the possibility to the future study.

4. Discussion

In this section, we discuss some important results we encountered in the previous sections. First, we want to study the meaning of the factors ${\alpha }_x$ and ${\alpha }_z$ appearing in Equation (9). A real dielectric material usually has ${\tilde{\epsilon }}_d>1$, so we get ${\alpha }_x=f+(1-f){\tilde{\epsilon }}_d>1$ and ${\alpha }_z^{-1}=f+(1-f)/{\tilde{\epsilon }}_d<1$, so both of them are larger than 1. Consider a harmonic $E$ field operating at a frequency $\omega$ much higher than ${\omega }_0$ and ${\omega }_p$. Then, according to Equation (17), the $Q$ field vanishes, and the $D$ field can be approximated as $[D_x,D_z]={\epsilon }_0[{\alpha }_x{E}_x,{\alpha }_z{E}_z]$. In one unit cell, the dielectric layer contributes $(1-f){\epsilon }_d{E}_x$ to the effective displacement field $D_x$, while the metal layer contributes $f{\epsilon }_0{E}_x$, so the sum of them gives $D_x=\left[f{\epsilon }_0+(1-f){\epsilon }_d\right]E_x={\alpha }_x{\epsilon }_0{E}_x$. This consideration also explains the form $Q_x=f{P}_x^m$, which is the dynamical part of the polarization field contributed by the metal layer.

To understand the origin of the ${\alpha }_z$ factor, we consider the situation of applying a time-varying electric field normal to the dielectric and metal layers. The applied $E$ field induces surface charges at the top and bottom surfaces of the dielectric layer, and the surface charges build an internal depolarization field inside the layer, pointing to the opposite direction of the applied field. This depolarization field cancels a part of the applied field, so we get an effective electric field smaller than the applied field. Using the constitutive relation $D_z={\epsilon }_0{E}_z+P_z={\epsilon }_0{E}_z+(1-f)P_z^d+f{P}_z^m$, we can derive the expression of $D_z$ in Equation (9) in a more straightforward manner:

$$D_z=\frac{D_z\left[D_z-(1-f)P_z^d\right]}{D_z-(1-f)P_z^d}=\frac{{\epsilon }_0{E}_z+f{P}_z^m}{1-(1-f)\left(1-\frac{1}{{\tilde{\epsilon }}_d}\right)}=\frac{{\epsilon }_0{E}_z+f{P}_z^m}{f+\frac{(1-f)}{{\tilde{\epsilon }}_d}}={\alpha }_z\left({\epsilon }_0{E}_z+f{P}_z^m\right).$$

(38)

This derivation also gives us the form of $Q_z={\alpha }_{z}f{P}_z^m$ automatically. This indicates that the depolarization effect also happens in the metal layer. In addition, the depolarization field plays the role of the restoring force acting to the “$Q_z$ oscillator” and gives a nonzero ${\omega }_0$ (see Equations (15), (23), (25) and (33b)).

The energy density can also be obtained by calculating the contributions from the dielectric layer and the metal layer separately and then summing them up. For example, the energy density corresponding to $E_x$ is given by

$$W_X=\frac{1}{2}\left(1-f\right){\epsilon }_d{E}_x^2+\frac{1}{2}f{\epsilon }_0{E}_x^2+fN\frac{m}{2}{\dot{x}}^2=\frac{{\epsilon }_0}{2}{\alpha }_x{E}_x^2+\frac{1}{2{\alpha }_x{\epsilon }_0{\omega }_p^2}{\dot{Q}}_x^2.$$

(39)

Here, $N$, $m$, and $x$ stand for the concentration, mass, and displacement of the conduction electrons in the metal, and the relations $P_x^m=Nqx$, $Q_x=f{P}_x^m$, and ${\omega }_{p0}^2=N{q}^2/m{\epsilon }_0$ have been used.

By the similar consideration, the energy density corresponding to $E_z$ can also be derived by considering the contributions from the dielectric and metal layers:

$$W_z=\frac{1}{2}\left(1-f\right){\epsilon }_d{E}_d^2+\frac{1}{2}f{\epsilon }_0{E}_z^{m2}+fN\frac{m}{2}{\dot{z}}^2=\frac{{\epsilon }_0}{2}{\alpha }_z{E}_z^2+\frac{1}{2{\alpha }_z{\epsilon }_{0}F{\omega }_0^2}\left({\dot{Q}}_z^2+{\omega }_0^2{Q}_z^2\right).$$

(40)

To obtain the right-hand side of the equation, we have used the relations: $P_z^m=Nqz$, $Q_z={\alpha }_{z}f{P}_z^m$, ${\omega }_{p0}^2=N{q}^2/m{\epsilon }_0$, and $E_z^m={\alpha }_z{E}_z+\frac{\left({\alpha }_{z}f-1\right)}{{\epsilon }_0}{P}_z^m$. The final relation is the consequence of the two different ways to write the displacement field: $D_z={\alpha }_z\left({\epsilon }_0{E}_z+f{P}_z^m\right)={\epsilon }_0{E}_z^m+P_z^m$.

When comparing the dielectric–metal multilayer system with the wire–SRRs array or the helical resonators array system [19,20], we find the most different feature of the present system is that it responds to the applied EM field anisotropically even at very high frequency, while the other two systems become isotropic in that limit. This difference stems from the fact that the relative permittivity 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_{\infty }+Q$ instead of $D={\epsilon }_{0}E+P$, and the prediction about the directional dependent behavior in the energy storage ability that was discussed in the previous section.

5. 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.