Causal homogenization of metamaterials

2.1 Deterministic causal homogenization of metamaterials

In the standard method of metamaterial homogenization (Figure 1) [3], we retrieve the effective refractive index n_eff and the impedance z_eff from the measured reflection coefficient r and the transmission coefficient t using the relations

Figure 1:

Homogenization of metamaterials.

(A) Schematic of a gold nanoparticle array. (B) The homogenized slab with the effective parameters. r and t are the reflection and transmission coefficients, respectively.

where the integer m is the branch, n₀ is the effective refractive index of the zeroth-order branch (m=0), k₀ is the wavenumber of the incident light, and d_eff is the effective thickness of the homogenized slab. In Eqs. (1)–(3), the undetermined integer branch m makes the effective parameters non-unique although the passivity condition requires that Re(z_eff)>0 and Im(n_eff)>0 [3]. Therefore, the homogenization procedure in metamaterial research admits multiple solutions for the homogenized permittivity and permeability.

In this work, we suggest that electromagnetic causality can determine the effective parameters uniquely. Causality is the most fundamental principle in physics because a physical action can occur only after its cause. In electromagnetics, the displacement field D(r, τ′) at the position r and time τ is determined only by the electric field E(r, τ′) for time τ′≤τ. For a linear medium, whose optical response is described by its permittivity ε(ω)/ε₀ relative to vacuum, the constitutive relation D(r, ω)=ε(ω)E(r, ω) determines the connection between D and E at the single frequency ω. If the constitutive relation of the fields D and E obeys causality, then the relative permittivity ε(ω)/ε₀ will have analyticity in the upper half ω plane [23]. This mathematical condition for ε(ω)/ε₀ gives rise to the KK relations

We may expect the same relation for the permeability μ(ω)/μ₀. However, it was noted that the KK relation applied to the permeability μ(ω)/μ₀, given by Eqs. (6) and (7), is not consistent with diamagnetism at zero frequency if the imaginary part Im{μ(ω)/μ₀} is positive for all ω [17]. Recently, it was argued that the imaginary part of permeability can be negative without violating physical laws [19]. Metamaterials with negative refractive index have effective parameters showing resonant and anti-resonant behaviors; that is, the imaginary part of permittivity is positive and that of permeability is negative near resonance [14]. On the other hand, we note that the induced magnetization is determined by the applied field B so that the causality should be applied to μ⁻¹. Nevertheless, the analytic properties of μ are also possessed by μ⁻¹ [17], i.e. μ⁻¹ is analytic in the upper half-plane if μ is analytic and has no zeros in the upper half-plane. Thus we adopt the same KK relation for the permeability μ(ω)/μ₀, namely

Here, we require the effective parameters of homogenized metamaterials to obey the causality and show that this requirement uniquely fixes the branch without ambiguity, thereby providing a deterministic parameter retrieval method. To impose the KK relation, we use the effective permeability μ_eff=n_effz_eff because the real and imaginary parts of n_eff and z_eff themselves do not obey the causality relation [24].

To obtain the causal condition for the branch m(ω) at the frequency ω, we first express the real and imaginary parts of the effective permeability in terms of the refractive index [Eq. (1)] and the impedance [Eq. (3)]. They are, respectively, given by

Plugging Eqs. (8) and (9) into the KK relation in Eq. (4), we obtain the Fredholm integral equation of the second kind for the integer branch m(ω) as follows:

where the functions g(ω) and K(ω′, ω) are given by

Equation (10) is the key result of this work. Previous works have applied the KK relation to the effective refractive index n_eff [8], [22]. However, the refractive index is not a physical response function respecting causality and there is no reason to require the KK relation [24]. Since Eq. (10) is a Fredholm integral equation of the second kind, it can be numerically solved. Discretizing the frequency by the N steps Δω′, Eq. (10) is converted to the matrix equation [25]

with the discretized frequencies ω_i and the identity matrix I of dimension N×N. By solving the matrix equation, Eq. (13), we can uniquely determine the branch integer m for the real part of the refractive index n_eff. Note that diagonal elements of the first matrix in Eq. (13) are zero because the principal values are defined to avoid the singularity of the integrand in Eq. (10).

The remaining ambiguity in determining the homogenized parameters is the effective thickness of the homogenized slab, d_eff. Metamaterials are usually composed of artificial nanostructures with uneven surfaces, and thus the boundary of the homogenized slab is often ambiguous. To determine the boundary of the homogenized slab, we assert that the effective thickness d_eff can be determined by minimizing the spectral averaged branch error

The bracket in Eq. (14) represents the spectral average over the domain of interest. The branch m should be integer in the ideal case, but homogenization artifacts such as the finite size of meta-atoms can cause non-integer solutions of m in the integral equation Eq. (10). Therefore, the averaged deviation of m from the integer in the spectral range of interest can be a measure of the homogenization artifacts. However, the following subsection will demonstrate that the choice of the effective thickness d_eff can minimize the deviation of m from an integer.

2.2 Homogenization of a gold nanoparticle array

To demonstrate the causal homogenization of metamaterials, we test the case of a gold nanoparticle array. Figure 2A shows a monolayer of gold nanoparticles of radius R=10 nm and periodicity P=26 nm, homogenized by the effective parameters. We vary the effective thickness d_eff to find the optimum value for minimizing the spectral averaged branch error Δm in Figure 2B. As the effective thickness d_eff increases, the branch m starts to deviate from the integer. The deviation is small at small frequencies but can be significant in the high-frequency region. Figure 2C shows that the spectral averaged branch error Δm is minimized at d_eff=2R=20 nm. Numerical simulation was carried out using the home-built finite-difference time-domain (FDTD) software. The grid size for the FDTD simulation is 1 nm. Two outer surfaces of the simulation domain normal to the wave propagation direction were perfectly matched layers (PMLs) of 10 grids, and four outer surfaces parallel to the propagation direction were determined by the periodic boundary condition (PBC). Material property of gold was taken from tabulated data [26]. The transmission coefficient t (the reflection coefficient r) is evaluated by the ratio between the incident fields and the transmitted (reflected) fields, which is spatially averaged at the plane normal to the incident wave vector in order to obtain the zeroth-order transmission (reflection). The wavelength range spans from 300 to 900 nm with the step size 1 nm (N=600), which corresponds to the frequency step size Δω′/2π=1.11 THz.

Figure 2:

Error minimizing choices of the effective thickness d_eff and the branch m.

(A) Schematic of a gold nanoparticle monolayer. (B) The branch m for different choices of effective thickness d_eff. (C) The spectral averaged branch error Δm against effective thickness d_eff.

Having determined the effective thickness, we can obtain the remaining effective parameters ε_eff, μ_eff, n_eff, and z_eff. Figure 3 shows the effective parameters of the homogenized gold nanoparticles. Near the frequency of 600 THz, we can find the surface plasmon resonance of the gold nanoparticles, which gives resonant peaks in the refractive index n_eff in Figure 3B. The permittivity and permeability can be simply obtained by the relations ε_eff=n_eff/z_eff and μ_eff=n_effz_eff, respectively. Note that the homogenized parameters presented in Figure 3 agree well with the theoretical predictions [27]. In Figure 3, our homogenization method for a monolayer yields m=0 because the monolayer of gold nanoparticles is too thin to have a nonzero branch m.

Figure 3:

Effective parameters for the homogenized gold nanoparticle monolayer.

(A) The real and (B) the imaginary part of the effective refractive index n_eff and (C) the real and (D) the imaginary part of the effective impedance z_eff of the monolayer (black lines) and seven layers (red lines) of gold nanoparticles.

To confirm the ability of our method to determine the nonzero branch m for bulky metamaterials, we homogenized seven layers of gold nanoparticles of periodicity P=26 nm in Figure 4. The homogenized slab is thick compared to the wavelength range of light. By minimizing the spectral averaged branch error Δm, we find that the effective thickness d_eff is 182 nm. Note that 7P=182 nm. As shown in Figure 4A,B, the effective impedance z_eff and the imaginary part of the refractive index Im(n_eff) are obtained uniquely because the homogenization procedure lacks ambiguity, as shown in Eqs. (1) and (3). The ambiguity occurs in the real part of the refractive index Re(n_eff). Figure 4C shows multiple branches of the real part of the refractive index Re(n_eff). Without considering the causality in the homogenization procedure, the branch m(ω) can be any integer that makes the refractive index continuous in the spectral window of interest. For example, the zeroth-order branch (m=0) is chosen for frequencies below 493 THz and the first-order branch (m=1) is chosen for higher frequencies to avoid abrupt, noncontinuous changes in the refractive index. Other choices are also possible with the integer difference Δm=1 at 493 THz. For example, we can choose m=1 below 493 THz and m=2 above 493 THz when the causality condition is not imposed. Another noncausal choice of m depicted in Figure 4C is m=−1 below 493 THz and m=0 above 493 THz. This degree of freedom when choosing the branch cannot be reduced without the causality condition. However, the causality condition can determine the set of branches m(ω) in the spectrum uniquely. Figure 4D shows the branch m solution for the integral equation Eq. (10) for seven layers of gold nanoparticles. The solution automatically guarantees the continuity of the real part of the refractive index Re(n_eff). We also confirm that the retrieved refractive indices of the seven layers are consistent with that of monolayer, as shown in Figure 3.

Figure 4:

Effective parameters for the homogenized seven layers of gold nanoparticles.

(A) The effective impedance z_eff and (B) the imaginary part of the refractive index n_eff of the seven layers of gold nanoparticles. (C) Examples of multiple branches of the real parts of the refractive index n_eff corresponding to the branch m without considering the causality condition. (D) The branch m determined by our homogenization method [Eq. (10)].