1. Introduction and background

Metamaterials are engineered composites, typically composed of multiple materials structured to enable unique properties distinct from those of their constituents [1]. Metamaterials have a variety of applications in cloaking, communications, sensing, microscopy, and optoelectronics [2,3]. To understand the properties of metamaterials, such as negative refractive index [3], and also because their structure and anisotropy renders their simulation computationally expensive, it is desirable to perform homogenization. This entails assigning a set of bulk properties that mimic the behavior of the structured metamaterial, as if it were a homogeneous slab [4–6].

Here, we explore the differences between various model classes that address metamaterials structured on wavelength and sub-wavelength scales, and propose a new homogenization method that can be applied to a variety of metamaterials. We suggest a robust solution to the branching problem in the popular Nicolson-Ross-Weir (NRW) homogenization method, which overcomes issues arising from the starting branch being determined under the assumption of electrical thinness. The method involves enforcing continuity after starting from the branch of the complex logarithm which has minimum absolute mean derivative at low frequencies. We also assign an optimal effective thickness to the slab. Subsequently, we compare the reflectance, transmittance, and absorptance of the original inhomogeneous material to the homogeneous slab that approximates it, using patterned films of PbS colloidal quantum dots as examples. Further, we discuss the Kramers-Kronig compliance of homogenized models, and observe that metamaterials with subwavelength inhomogeneities are Kramers-Kronig compliant but those with larger inhomogeneiteies may not be. In this work, we define “effective” and “equivalent” as they are used in [7,8]. “Effective” refers to a physical homogenization in which the wavelength is much greater than the feature size of the inhomogeneities, and, as a result, the field within the effective slab will be similar to the macroscopic field in the heterogeneous medium. “Equivalent” is a less restrictive term, which only attempts to reduce the metamaterial to a homogeneous slab with the same reflectance, transmittance, and absorptance. If the size of the inhomogeneities is comparable to the wavelength, an effective medium is not achievable; however an adequate equivalent medium can still be found [7]. An example application of this method is in the incorporation of a homogeneous model into the simulation of a larger composite device structure, allowing for the use of a fast 1D simulation framework such as the Transfer Matrix Method rather than a costly 2D-3D full-wave simulation. A demonstration of this is provided in [9].

1.1 Effective and equivalent materials parameter extraction models

There are a variety of models for assigning effective or equivalent parameters to materials. A simple option is to take a volume average of the component refractive indices. The Maxwell Garnett mixing formula [10] uses a point-dipole approximation to obtain an expression for the permittivity for arbitrarily shaped particles in a medium, provided their volume fraction is small. This constraint limits the applicability to a narrow range of systems [10]. The Bruggeman mixing model provides a more symmetric treatment of host medium and inclusions, and is thus applicable for any volume fraction. However, the Bruggeman model also does not account for geometry [10]. Likewise, considering particles as points without regard for particle geometry can create models that diverge [11] or violate basic principles such as causality [12]. These methods also often break down in the presence of strong resonances [11]. Both methods assume inhomogeneities large enough to be described by a macroscopic permittivity but still smaller than the wavelength [13], making them effective medium models. Other approaches have also been proposed [8], such as using average field strengths [14], Drude-Lorentz models [15], or Bloch mode analysis [16].

Among the most popular homogenization methods is the Nicolson-Ross-Weir (NRW) method [17,18]. This involves obtaining measurements of $S_{11}$ and $S_{21}$, which are the backward (reflected) and forward (transmitted) complex scattering coefficients, respectively. This can be done experimentally or via computational simulation. The original paper by Nicolson and Ross uses a TEM-mode fixture to obtain time-domain measurements of the S-parameters at microwave frequencies [17], and Weir proposes an analogous method that directly takes frequency-domain measurements [18]. For optical regimes, to obtain the S-parameters, various interferometric methods (e.g. Michelson interferometry) can determine the phase difference between the incident and reflected or transmitted beams, thus allowing determination of the complex S-parameters [19–21]. In practice, if there is a substrate, it should also be accounted for. The values of $S_{11}$ and $S_{21}$ are then used to obtain the complex refractive index and wave impedance, from which permittivity and permeability can be obtained. While the imaginary part of the refractive index can be determined exactly, the equation for the real part of the refractive index involves a complex logarithm with multiple branches, thus resulting in a branching problem. The modal analysis, Maxwell Garnett, and NRW methods are compared in [22,23] and their agreement is used to justify the homogenization procedure.

1.2 Nicolson-Ross-Weir method

In the Nicolson-Ross-Weir method [4,24,25], the S-parameters are given as:

The sign in the equation for $Z_{\text {eff}}(\omega )$ is chosen to ensure that $\text {Re}\{Z_{\text {eff}}(\omega )\}$ and $\text {Im}\{N_{\text {eff}}(\omega )\}$ are both nonnegative, which is the same as enforcing $|e^{iN_{\text {eff}}(\omega )k_0d_{\text {eff}}}| \leq 1$ due to passivity considerations [4,24]. Equation (3) can be rewritten as:

Thus, knowledge of the refractive index and wave impedance is equivalent to knowledge of the permittivity and permeability.

1.3 Past solutions to the branching problem

Many approaches have been taken to solve the branching problem. One approach is via the Kramers-Kronig relations, given by:

A different approach is to use $\text {d}n_\text {eff}(\omega )/\text {d}\omega$ to enforce continuity of $n_\text {eff}(\omega )$ [29], or to use its Taylor expansions and take the initial branch consistent with passivity [24]. A more rigorous treatment involves the analytic continuation of the logarithm [25,30,31]. An alternative method for choosing the initial branch number extends the starting frequency for a small number of metamaterial layers to a larger number [5]. Similarly, assuming electrical thinness at a sufficiently low frequency results in a small exponent in Eq. (3) and thus $m = 0$ being the starting branch [32]. However, this requires the wavelength in the sample to be known, which is often not the case, to ensure frequencies below the first branch transition are included [7]. If an estimate is used, the method will lack robustness, an issue we address with this work. For thick slabs, one can compare two thicknesses to find the branch on which they agree; however, the problem of needing to know the wavelength in the sample persists [7]. Finally, an alternative solution to the branching problem, using deep learning, is presented in [33].

1.4 Consideration of the Kramers-Kronig relations

The NRW method has been reported to yield materials parameters violating basic passivity or causality conditions [8]. Moreover, in the presence of nonlinear and saturated polarizability, Eq. (7a) may again fail to hold [13,34]. It is also suggested that in systems with gain, the refractive index is not required to satisfy a Kramers-Kronig relation (but its square is) [34]. When the averaged electric fields are the same for the original material and its effective counterpart (i.e. when the wavelength is much greater than the scale of the inhomogeneities) the Kramers-Kronig relations must hold as a result of causality [12]. However, in large-feature regimes, where the only considerations are the macroscopic reflectance, transmittance, and absorptance spectra of the original metamaterial and its equivalent slab, this is no longer a requirement, as demonstrated in the present work. The latter approach to homogenization is what we refer to as an “equivalent” model rather than an “effective” one, consistent with [7] and [8]. This does not, however, suggest that causality is violated, or that superluminal information transfer is possible [12]. Rather, the equivalent model is not physical: it is simply the model that results in the same output reflectance, transmittance, and absorptance as the original inhomogeneous material. It is unphysical to treat its internal workings as homogeneous since the wavelength is on par with the feature size and effects such as diffraction and scattering dominate. A related discussion of causality-violating models produced when using point-dipole approximations is available in [12]. In these cases, the true geometry approximated by a point dipole interacts with a wave before the imagined dipole at its center, thus making it seem as if the response to the wave occurred before the interaction with the object.

1.5 Effective thickness

The effective thickness of a metamaterial is distinct from its geometric thickness. In [24], the effective thickness is calculated via the boundaries at which incident and outgoing waves are planar. The effective thickness can also be handled by rounding the branch number obtained from the Kramers-Kronig relations and taking the effective thickness as the number that minimizes the rounding error [6]. We may make a further argument regarding how interaction with geometric structures prior to reaching an assumed point dipole creates models that violate causality [12]. Since the plane-wave behavior breaks down significantly for equivalent materials, it is possible that this interaction before the geometric boundary is responsible for the violation of the Kramers-Kronig relations. This is discussed further in Section 3.3.

1.6 Motivation

Here, we construct a method to create effective and equivalent models for inhomogeneous materials that solves the branching problem in the NRW method. The process involves a simple approach of starting at the flattest branch and employing continuity. We use the assumption that any medium will act as a homogeneous slab at very long wavelengths and thus eventually the materials parameters will stabilize to the DC limit.

Together with the consideration that the branches will separate for low frequencies, the unphysical branches will steepen while the correct one will not. Even if the starting branch index is not zero, we can use continuity to trace through the rest of the spectrum. This notion of continuity is that suggested in [7] and [28]. Thus, the “hard” requirement that the sample be thin enough to start at a branch index of zero is eased to a “soft” requirement that the correct branch have a lower absolute mean derivative than the others, a simple and intuitive requirement true for a wide range of materials. It is in this sense that our method is robust, requiring only a comparison across a range of branches to find the one with minimum absolute mean derivative. We confirm the validity of our results via the analytically calculated optical properties, observing high agreement with the FDTD simulation of the inhomogeneous slab. Moreover, we also minimize the error in these comparisons to determine an effective thickness for the slabs. We argue that effective models obey the Kramers-Kronig relations since it is physically justifiable to treat them as homogeneous while equivalent models need not obey the Kramers-Kronig relations due to the models’ unphysical nature, where they are simply used to match optical behavior.

2. Methods

2.1 Solution to the branching problem

We use the Lumerical Finite Difference Time Domain (FDTD) solver [35] to simulate and extract the scattering coefficients $S_{11}$ and $S_{21}$ for the inhomogeneous material, and subsequently make all calculations up to and including Eq. (5a) and Eq. (5b).

Our resolution is to start from the branch with the minimum derivative in the low frequency limit, then use continuity to trace the correct function through the branches. At sufficiently low frequencies, optical effects associated with finite inhomogeneities appproach their low frequency limit behavior. Here, the refractive index, permittivity, and permeability are constants for a wide variety of materials.

We start by using approximate physical considerations to set an upper branch limit. Branch number is connected to the number of complete periods traversed by the wave in the material. We are mainly concerned with the visible and near infrared (NIR) portions of the spectrum, so we consider a rough upper limit of 10 for the refractive index and a lower wavelength limit of 200 nm, leaving some tolerance. Then, to obtain the approximate number of complete periods traversed in a round trip in the material, we can write:

The branch with the lowest absolute mean derivative can be chosen since the term in Eq. (5a) containing the branch index is a reciprocal of $k_0$ and thus of the frequency, meaning that the false branches become very steep for low frequency (to be more precise, for a small value of $k_0 d_{\text {eff}}$, which is dimensionless). Then, in order to find the minimum absolute mean derivative at low frequencies, we may denote the $m^\text {th}$ branch of $n_\text {eff}$ as $n_{\text {eff},m}(\omega _j)$, where $j$ is the index denoting the $j^\text {th}$ frequency point considered, starting from the lowest ($j = 1$). Then for a suitably chosen $j'$ such that a discontinuity is avoided, we may write the branch index for the first point $m_1$ as:

Subsequently, we may retrieve the remaining branch indices via the assumption of continuity, using the index from the current branch and its two nearest neighbors and selecting the closest one. We have also attempted a first-order Taylor series based method for this, but it did not produce a tangible difference and is thus not presented in this work. The branch index $m_{j+1}$ is found from $m_{j}$ via the recursive relation:

2.2 Effective thickness and modeling the effective/equivalent slab

Next, we calculate the reflectance $R$, transmittance $T$, and absorptance $A$ of both the original heterogeneous material and its effective or equivalent model. The inhomogeneous values $R_\text {inh}(\omega )$, $T_\text {inh}(\omega )$, and $A_\text {inh}(\omega )$ were calculated via FDTD directly, using a plane-wave source. The corresponding values for the effective and equivalent media were calculated via the theoretical consideration of a single slab in free space. The solution for the optical behavior of a single slab is well-known and available in [36]. We may write:

These are essentially the same formulae that we used when calculating the effective parameters through the NRW method and the S-parameters. Thus, we expect agreement to be exact in the case of an effective medium, and approximate in the case of an equivalent medium, which will fail to account for all the relevant optical effects.

We further refine our method to find the effective thickness in the wavelength range of interest, taking advantage of the slab equations presented above. We optimize for the value of $d_{\text {eff}}$ by using the NRW method and the associated calculation of the reflectance, transmittance, and absorptance and calculating the mean squared error (MSE) between our FDTD simulations and slab calculations as below:

 

Fig. 1. Schematic of the procedure employed in this work for metamaterial homogenization, via optimizing the NRW method over the effective thickness. Following a single FDTD simulation to obtain the S-parameters, the modified NRW method is applied for different candidates for the effective thickness ($d_0$, $d_1$, …, $d_\text {max}$), and the $d$ with the lowest mean squared error (MSE) compared to the FDTD calculation is chosen.

Download Full Size | PDF

3. Results and discussion

We demonstrate the accuracy of our method with a variety of examples. We first check that our method is self-consistent and accurately retrieves parameters from a homogeneous slab. We then move to patterned thin films of PbS colloidal quantum dots (CQDs), and demonstrate that the equivalent models we create are accurate, and that this accuracy is further improved with effective thickness optimization. Afterwards, we compare two different size regimes of the same structure, one corresponding to an equivalent model and the other to an effective model (with respect to the visible and NIR wavelengths). We then investigate the Kramers-Kronig compliance of these models. Finally, we introduce an example of a birefringent grating structure and compare the spectra obtained at different polarizations. The materials we use have applications to hierarchically structured spectrally selective optoelectronics [37].

3.1 Homogeneous slab example

As an example, we start by considering a 790 nm thick homogeneous slab of Si. The calculated branches, the branch number, the extracted refractive index, and the optical behavior are shown in Fig. 2. For homogeneous slabs, the effective thickness and geometric thickness are unsurprisingly the same. The agreement in optical behavior is excellent, and any error is explained by FDTD simulation resolution. The resulting MSE is $4.69 \times 10^{-5}.$ The NRW formulas account for interference effects due to interfaces, and thus the extracted refractive index is independent of these effects, as can be verified from this example of a homogeneous slab. The extracted solution is the correct combination of the branches, as can be seen in Fig. 2. Most notably, the first branch index is $m_1 = 1,$ demonstrating the capability of the method to start from a nonzero branch index and thus removing the previous strict requirement of a thin slab [32]). The ability of our model to correctly start from nonzero branch numbers demonstrates its robustness.

 

Fig. 2. Demonstration of the method using the example of a homogeneous slab of Si. (a) Branches of the complex logarithm (dotted) and real ($n_\text {eff}$) and imaginary ($\kappa _\text {eff}$) parts of the extracted refractive index (solid). (b) Variation of branch number with frequency. (c) Extracted refractive index (solid lines) and input refractive index (dashed lines) of the Si slab. (d) Reflectance (blue), transmittance (green), and absorptance (red) from the FDTD simulation of the inhomogeneous slab ($R_\text {inh}, T_\text {inh}, A_\text {inh}$; solid lines) and the analytic homogeneous slab formulas ($R_\text {eff}, T_\text {eff}, A_\text {eff}$; dashed lines), showing agreement.

Download Full Size | PDF

3.2 Equivalent models for patterned slabs

We next demonstrate the model for an array of PbS CQD disks in free space, a patterned structure relevant in spectrally selective optical devices such as solar cells [37]. The disk radius is 253 nm, the period is 632 nm, and the thickness is 790 nm. PbS CQDs were synthesized using a modification of the Hines Method [38,39], with a nominal target exciton peak of 1240 nm. Samples were prepared by spin-casting the PbS CQD solution in octane onto a roughly 25 mm $\times$ 25 mm bare silicon wafer. Approximately 40 $\mu$l of CQD solution was drop-cast (through a 0.22 $\mu$m filter) onto the bare silicon wafer and spun at 2000 RPM for 30 seconds. This produced a roughly 100 nm layer of oleic acid-capped PbS CQDs. Variable Angle Spectroscopic Ellipsometry (VASE) measurements were taken using a J. A. Woollam VASE UV-NIR system over a wavelength range of 300 to 2500 nm. Ellipsometry constants were fit using a Kramers-Kronig consistent general oscillator (GenOsc) model with the WVASE software [40]. The optical model and thickness converged at a film thickness of 91 nm, agreeing well with the expected value from the procedure above. The complex refractive index functions extracted by this method can be seen in the inset of Fig. 3(b).

 

Fig. 3. The use of the method for a patterned slab. (a) The inhomogeneous structure, a PbS CQD thin film patterned into an array of disks. The disk radius is 253 nm, the period is 632 nm, and the thickness is 790 nm. (b) The equivalent refractive index model, with real ($n_\text {eff}$) and imaginary ($\kappa _\text {eff})$ parts shown. Inset shows the real (dashed blue line) and imaginary (dashed red line) parts of the refractive index of a homogeneous slab of PbS CQDs obtained via VASE measurements. (c) The MSE as a function of effective thickness, with geometric thickness of 790 nm and optimal effective thickness of 948 nm. (d) Reflectance (blue), transmittance (green), and absorptance (red) from the FDTD simulation of the inhomogeneous slab ($R_\text {inh}, T_\text {inh}, A_\text {inh}$; solid lines) and the analytic homogeneous slab formulas ($R_\text {eff}, T_\text {eff}, A_\text {eff}$; dashed lines).

Download Full Size | PDF

The results are given in Fig. 3. Note that since the inhomogeneities are of size comparable to the wavelength, this is only an equivalent model for most of the wavelength range of interest. This explains the mild discrepancies between the inhomogeneous structure and its equivalent model, since the latter cannot account for effects such as diffraction. It is worth noting that the peak wavelength locations mostly match between the two sets of spectra in Fig. 3(d); however, the amplitudes are slightly different. For sake of comparison, Fig. 4 displays the optical behavior with and without effective thickness optimization. As can be seen in the figure, there is a value of $d_{\text {eff}}$ minimizing the MSE that is slightly thicker than the slab itself, supporting the notion that an extra buffer is needed around the material to account for non-plane wave behavior. We thus demonstrate the improvement in our model that effective thickness optimization provides.

 

Fig. 4. Comparison of the reflectance (blue), transmittance (green), and absorptance (red) from the FDTD simulation of the inhomogeneous slab ($R_\text {inh}, T_\text {inh}, A_\text {inh}$; solid lines) and the analytic homogeneous slab formulas ($R_\text {eff}, T_\text {eff}, A_\text {eff}$; dashed lines) under the (a) presence and (b) absence of effective thickness optimization for the patterned slab in Fig. 3.

Download Full Size | PDF

3.3 Kramers-Kronig compliance of effective and equivalent models for patterned slabs

We next demonstrate two size regimes of the aforementioned patterned PbS CQD structure to compare an effective medium (with feature size much smaller than the wavelength) with an equivalent one (with feature size comparable to the wavelength). The precise geometry of the former consists of a 158 nm period and 63 nm disk radius, with a thickness of 790 nm, while the latter is the example of Fig. 3. The refractive index models and spectra for the two regimes are shown in Fig. 5. We compare $n_{\text {eff}}$ for each structure to that found by applying the Kramers-Kronig relation of Eq. (7a) to $\kappa _{\text {eff}}$. The results are also given in Fig. 5. Slight discrepancies in the Kramers-Kronig relations are expected due to truncation of the integral and the numerical handling of the point in the integral that would normally cause divergence. We observe that, for the effective model, the Kramers-Kronig relations are obeyed within bounds of the aforementioned error. However, there is major discrepancy with the predictions of the Kramers-Kronig relations in the case of the equivalent model. This is explained by the unphysical nature of the homogeneity approximation on this scale. Furthermore, the real part of the refractive index calculated from the Kramers-Kronig relations and that calculated from the NRW method approach each other for longer wavelengths, where feature size becomes effectively smaller in comparison, and the model moves towards the “effective” regime. This behavior further supports the connection between Kramers-Kronig compliance and the effective regime. We also note that the agreement between the modeled and simulated slabs is almost exact for the case of feature size much smaller than the wavelength, but less so in the case where wavelength and feature size are comparable. This discrepancy is alleviated by the optimization for effective thickness, as discussed previously.

 

Fig. 5. Comparison of effective (158 nm period and 63 nm disk radius) and equivalent (632 nm period and 253 nm disk radius) PbS CQD patterned slab arrays (both thicknesses are 790 nm). For the material with feature size much less than wavelength, the refractive index model (imaginary part $\kappa _\text {eff}$, real part via Kramers-Kronig transformation $n_\text {eff, KK}$, and real part via the algorithm in this work $n_\text {eff}$) is given in (a) and comparison of the reflectance (blue), transmittance (green), and absorptance (red) from the FDTD simulation of the inhomogeneous slab ($R_\text {inh}, T_\text {inh}, A_\text {inh}$; solid lines) and the analytic homogeneous slab formulas ($R_\text {eff}, T_\text {eff}, A_\text {eff}$; dashed lines) is given in (b). Similarly, for the material with feature size comparable to wavelength, refractive index model and optical behavior are given in (c) and (d), which are repeated from Fig. 3 for ease of comparison.

Download Full Size | PDF

Finally, we compare the field profiles for the two cases (Fig. 6), to illustrate the plane wave behavior and how it relates to the effective thickness and adherence to the Kramers-Kronig relations. As can be seen, plane wave behavior persists until almost the geometric boundary of the material for the case of the effective medium model, justifying the effective thickness being close to the geometric thickness for effective models; however, for the equivalent case, the plane wave behavior begins and ends further from the boundaries. As visible in Fig. 3, the effective thickness is 948 nm, and it can be seen in Fig. 6 that the non-plane wave behavior is indeed approximately within such a range, rather than the geometric thickness of 790 nm. Additionally, we note that the effective thickness is calculated to account for the behavior across all frequencies, so it includes a consideration of where plane wave behavior starts for each frequency, not just those shown.

 

Fig. 6. Comparison of the magnitude of the electric field calculated using FDTD simulations at 1060 nm wavelength for (a) the PbS CQD disk array with feature size much less than (158 nm period and 63 nm disk radius) and (b) comparable to (632 nm period and 253 nm disk radius) the wavelength (both thicknesses are 790 nm). In both cases, the plane wave excitation source, with amplitude $E_0$, is incident from the bottom of the image.

Download Full Size | PDF

3.4 Birefringent grating

We also demonstrate the method on an anisotropic grating, with differences in the refractive index spectra at parallel and perpendicular polarizations of the electric field relative to the grating. The grating, made of strips of a PbS CQD thin film, is infinite in one direction, with a width of 200 nm, thickness of 790 nm, and a periodicity of 400 nm. The results are given in Fig. 7. The models are both equivalent models due to the feature size, and as expected we observe slight disagreement in optical spectra and a lack of Kramers-Kronig compliance. The birefringence, defined as $|N_{\text {eff},\parallel } - N_{\text {eff},\perp }|$, is found as 0.35 at 2000 nm with $N_{\text {eff},\parallel } = 1.62$ and $N_{\text {eff},\perp } = 1.27$. It is worth noting that simpler models that do not account for geometry (such as volume averaging) would not be able to make such a distinction in refractive index for birefringence.

 

Fig. 7. Demonstration of the method proposed in this work on a birefringent grating, which is infinite in one direction, with a width of 200 nm, thickness of 790 nm, and a periodicity of 400 nm. With electric field polarization (blue arrow) parallel to the grating (a), the equivalent refractive index (imaginary part $\kappa _\text {eff}$, real part via Kramers-Kronig transformation $n_\text {eff, KK}$, and real part via the algorithm in this work $n_\text {eff}$) is shown in (c) and comparison of reflectance (blue), transmittance (green), and absorptance (red) from the FDTD simulation of the inhomogeneous slab ($R_\text {inh}, T_\text {inh}, A_\text {inh}$; solid lines) and the analytic homogeneous slab formulas ($R_\text {eff}, T_\text {eff}, A_\text {eff}$; dashed lines) is shown in (e). Under perpendicular polarization (b), the equivalent refractive index is shown in (d) and the corresponding optical spectra in (f). The blue arrow represents the electric field, and the green arrow represents the magnetic field. The incident wave propagates from behind the structure and out of the figure.

Download Full Size | PDF

4. Conclusion

In this work, we have proposed and demonstrated a method for metamaterial homogenization and parameter extraction that robustly solves the branching problem in the NRW method and incorporates effective thickness. We also analyze our results in the context of Kramers-Kronig relations and elucidate the difference between effective and equivalent models.

Our method introduces the idea of starting at the branch of the complex NRW logarithm which begins with the minimum absolute mean derivative, and then enforcing continuity of the function to obtain the rest of the real part of the refractive index. This method is robust to the frequency range used compared to previous models. Moreover, we compare the reflectance, transmittance, and absorptance of the heterogeneous model to that of the effective or equivalent homogeneous slab and observe nearly exact agreement for the former. This indicates that the homogenization process is physically valid since feature size is much smaller than wavelength. There are mild discrepancies for equivalent models where effects such as diffraction cannot be accounted for by a homogeneous slab, but the optical behavior is nonetheless close enough for most practical purposes. We choose the effective thickness of the modeled slab such that the deviation from the optical behavior of the original material is minimized. The effective thickness is the same as the geometric thickness for an effective model but typically thicker for an equivalent model due to the extended region in which the fields do not behave as a plane wave. Thus, our comparison of optical behavior not only serves as verification of our model, but also allows improvement of the results by enabling us to choose the optimal effective thickness. We have also examined the compliance of the effective and equivalent medium models with the Kramers-Kronig relations. We show that effective models are Kramers-Kronig compliant due to the physical nature of assuming homogeneity, whereas equivalent models do not necessarily comply with the Kramers-Kronig relations due to the somewhat unphysical nature of the homogeneity assumption, and not due to an actual violation of causality. This noncompliance does not invalidate the model, since a homogeneous model with the same optical behavior is still produced. Future work demonstrating that this method can also be used for the extraction of anisotropic permittivity tensors and comparing the results to experimental measurements of example metamaterials is also planned.

In conclusion, our compact method can obtain both effective and equivalent material models, which greatly simplifies both simulation and interpretation of heterogeneous media in complex photonic structures. This method could be deployed across a variety of applications that integrate multiple components and are computationally expensive to simulate since it allows for solving or optimization by 1D means such as the Transfer Matrix Method rather than 2D or 3D methods such as full-wave simulations [9], broadening the applications for metamaterial components.