1. Introduction

Over the recent years, a constantly escalating interest has been perceived for metamaterials, a class of artificial media which exhibit negative constitutive parameters at certain frequency bands. A typical realization of these substances lies on periodic arrays of thin metallic rods and split-ring resonators (SRRs) with dimensions much smaller than the operating wavelength [1]-[4], or properly interconnected lumped circuit particles [5], [6]. Up till now, a great deal of wave guiding and radiation components that exploit the unusual, yet challenging, properties of metamaterials (e.g. antiparallel phase and group velocities or negative refraction) has been proposed, for advanced bandwidth features and smaller packaging sizes [7]-[14]. Along with them, one should mention the impressive superlensing effect and the amplification of evanescent waves [15]-[21], as well as the wire medium lenses, that achieve superior resolution via the canalization of electromagnetic waves through them [22], [23]. Actually, most of these characteristics owe their existence to the interactions between regular media and metamaterials, which give rise to resonant surface polaritons traveling toward parallel directions to the interfaces.

In essence, surface polaritons are modes appearing at the interfaces between double positive (DPS) and double negative (DNG) or single negative (SNG) materials, which undergo an exponential decay on both sides of the interface [24]. On the other hand, metamaterial slabs, also, support conventional surface waves, with a standing wave profile, in their interior [25], [26]. Taking advantage of their resonant nature and their property to focus the electromagnetic energy around the interface, an assortment of microwave devices, like waveguides and patch antennas with SNG/DNG inclusions and dimensions far below the diffraction limit, have been developed [7], [11]. Furthermore, the form of the surface waves dispersion diagram in grounded metamaterial slabs, that results in bandgap regions where no surface wave is supported, motivated the use of these structures as patch antenna substrates to subdue parasitic radiation and accomplish enhanced directivity [8]. However, all previous contributions employ an effective metamaterial that presumes homogeneous constitutive parameters, thus not examining the significant influence of the finite-sized metallic inclusions that construct practical DNG media. To this end, an experimental study dealing with the magnetic surface polaritons of SRR-based materials through an attenuated total reflection technique has proven that such modes are, indeed, excited [27]. Indirect evidence, also, for the evolution of surface polaritons can be found in [17], [28], where evanescent wave amplification is computationally and experimentally explored.

It is the purpose of this paper to present a novel analytical and numerical study of surface waves, including resonant surface polaritons, that emerge at various types of realistic metamaterial configurations. More specifically, it is investigated in what degree the local homogenized constitutive parameters of various metamaterial arrangements can be used for the description of the surface waves characteristics of real non-homogenized structures. Regarding the analytical framework, complete dispersion diagrams are extracted by incorporating the average parameters of planar SNG and DNG slabs - engineered by periodical patterns of broadside SRRs and S-rings - in the corresponding dispersion relations. Despite the existence of fully analytical formulae for the average electric permittivity and magnetic permeability, their evaluation via a set of periodically-realized finite difference time domain (FDTD) simulations is preferred, to drastically improve the overall accuracy. In this context, the non-local effective constitutive parameters are, first, computed through the unit cell’s reflection and transmission coefficients and, then, the average parameters are obtained via the technique derived in [29]-[34]. Successively, surface waves dispersion characteristics are investigated on an elaborate basis by means of an efficient numerical algorithm, which involves the precise FDTD modeling of the preceding structures and the application of consistent periodic boundary conditions at the domain edges. The results, so acquired, are compared to those of the analytical scheme and instructive deductions about the important role of the DNG constitutive elements finite size are systematically drawn. It is emphasized, herein, that the principal difference between the analytical and the numerical formulation is that the former assumes a local response for the artificial medium, described by its average parameters, which, nonetheless, may become inaccurate as the wavevector approaches the limits of the first Brillouin zone.

2. Development of the theoretical methodology

The key issue of the proposed formulation is based on a rigorous effective medium approach for the elaborate investigation of surface modes attributes, evolving at the interfaces of planar DNG/SNG arrangements. In particular, the necessary surface wave dispersion diagrams for structures composed of thin wires and SRRs are extracted through the macroscopic rendition of constitutive parameters, defined via the respective electric and magnetic field intensity mean values. The geometry under examination is depicted in Fig. 1, with d the width of the slab along z-axis and ε̿(ω), μ̿(ω) the electric and magnetic permeability tensors for the description of the metamaterial. Following a precise electromagnetic analysis, it can be proven that such slabs may guide TM_x (magnetic field intensity perpendicular to x-axis) and TE_x (electric field intensity perpendicular to x-axis) polarized waves, parallel to their faces, whose dispersion relations between frequency ω and the y-oriented wavenumber component k_y are given by

 

Fig. 1. An infinite planar DNG slab with effective constitutive parameters ε̿(ω) and μ̿(ω).

Download Full Size | PDF

In the above, k′z = [(ω/c)² μ_yy ε_xx − (μ_yy/μ_zz)k_y ²]^1/2 and k′z = [(ω/c)² μ_xx ε_yy − (ε_yy/ε_zz)k_y ²]^1/2 are the z-directed wavenumber components inside the slab for the TM_x and TE_x case, respectively, a_z = (k_y ² − ω ²/c ²)^1/2 denotes the z-directed attenuation factor in the air, c is the speed of light in the air, while the plus and minus signs indicate symmetrical and antisymmetrical modes with regard to the z = 0 plane. Bear in mind that the propagation constant along x-axis has been assumed zero, since for k_x ≠ 0, the wire medium, considered parallel to x-axis, exhibits nonlocal effects, which are beyond the scope of the present investigation. Additionally, both constitutive parameter tensors are selected to be of the ε̿(ω) = diag[ε_xx,ε_yy,ε_zz] and μ̿(ω) = diag[μ_xx,μ_yy,μ_zz] diagonal form, as in the thin-wire/SRR structure. In this context, (1) and (2) reveal that surface polaritons, with an exponentially decaying profile in the interior of the slab, appear, solely, at frequencies where μ_yy < 0 (TM_x case) or ε_yy < 0 (TE_x case).

Prior to the solution of (1), (2) and the extraction of the dispersion characteristics for the supported surface waves, it is necessary to determine the type of the ε̿(ω) and μ̿(ω) components. Starting from wire media with a parallel orientation toward x-axis and presuming that propagation occurs on the yz-plane, one can reliably describe the xx-component of the electric permittivity tensor by means of a Drude model relation, as

 

Fig. 2. Geometry of (a) an edge-coupled SRR, (b) a broadside SRR, and (c) an S-ring.

Download Full Size | PDF

where the plasma frequency ω_p depends on the wire and the unit-cell dimensions. On the other hand, the edge-coupled SRR, shown in Fig. 2(a), accomplishes a magnetic activity along y-axis, which complies with the Lorentz-type relation

for a a constant and ω _m0 the resonance frequency, both specified by the geometry of the unit cell. Moreover, such an SRR displays an electric activity toward z-axis with a constant effective electric permittivity ε_zz. Similar behavior is observed for the broadside SRR and the S-ring, presented in Figs 2(b) and 2(c), respectively, except for the fact that ε_xx in the latter setup follows a Lorentz-type model. Taking into account the above remarks, the slab must contain SRRs on constant z-planes to attain μ_yy < 0 and thus allow the excitation of TM_x polarized surface polaritons. In contrast, TE_x polaritons can not evolve because neither the SRR structure nor the x-oriented wire medium result in negative ε_yy.

Although analytical expressions for ω_p, ω _m0, a, appearing in (3) and (4), may be found in the relevant literature [1]-[4], we, herein, prefer to calculate them via proficient FDTD simulations and therefore significantly enhance the levels of total accuracy. The main reason for this choice is that existing formulae are, only, approximate, since they are acquired by means of simplified assumptions for the current and charge distribution on the surface of the metallic metamaterial particles. The process of deriving the necessary average constitutive parameters involves, as its first step, the calculation of effective parameters ε _eff and μ _eff for an one-cell slab (i.e. a slab infinite along the y-axis and one-cell thick along the z-axis), through the reflection and transmission coefficient for normal incidence [29], [30]. Notice that the parameters, so extracted, are non-local and valid solely for normal incidence. Moreover, if electric and magnetic resonances exist at close frequencies - as in the case of the S-ring arrangement - ε _eff and μ _eff are applicable only to the slab width for which they have been obtained [33], [34]. Subsequently, following [32], the required average parameters ε̅ and μ̅ can be evaluated in terms of

with θ the phase advance along the unit cell and S_b a parameter equal to 1 or -1 depending on the electric or magnetic resonances that dominate the metamaterial.

For the FDTD analysis to be precise, periodic boundary conditions are applied along the x-and y-axis in an effort to model the infinite dimensions of the slab parallel to the corresponding directions. Essentially, these conditions are of the ψ(r+p) = ψ(r)exp(−j q·p) form, where ψ is any field component, r denotes the position vector, p = p_x x̂ + p_y ŷ the primary lattice vector and q = q_x x̂ + q_y ŷ the Bloch vector. In fact, since for the extraction of the effective parameters our interest focuses on a normally impinging plane wave, it is deduced that q_x = q_y = 0. Furthermore, p_x, p_y are the dimensions of the unit cell along the x, y axes, respectively.

 

Fig. 3. Numerically extracted average constitutive parameters ε̅_xx and μ̅_yy of (a) the broadside SRR with w = 3 mm, s = 0.25 mm, g = 0.5 mm, h = 0.5 mm, p_x = p_y = 5 mm, L = 5 mm and (b) the S-ring with w = 0.5 mm, s = 0.25 mm, h = 0.5 mm, ℓ₁ = ℓ₂ = 3 mm, p_x = 5.5 mm, p_y = 2.5 mm, L = 4 mm.

Download Full Size | PDF

In this context, Fig. 3(a) illustrates the average permittivity ε̅_xx and the average permeability μ̅_yy of the broadside SRR (Fig. 2(b)) with w = 3 mm, s = 0.25 mm, g = 0.5 mm, h = 0.5 mm, p_x = p_y = 5 mm, and L = 5 mm, where L denotes the unit cell length along the z-axis. To compute ε̅_xx and μ̅_yy, we employ (5) and (6) with S_b = -1, because of the magnetic SRR resonant behavior. From the plots, one may, readily, discern the Lorentzian resonant behavior of μ̅_yy around 11.68 GHz and its negative values between 11.68 GHz and 12.37 GHz. Conversely, ε̅_xx receives a continuous, yet slow, augmentation with frequency, thus allowing its treatment as a constant for narrow bands. A similar performance can be, also, acquired for the edge-coupled SRR (Fig. 2(a)). Concerning the S-ring (Fig. 2(c)), its average parameters are shown in Fig. 3(b), for w = 0.5 mm, s = 0.25 mm, h = 0.5 mm, ℓ ₁ = ℓ ₂ = 3 mm, p_x = 5.5 mm, p_y = 2.5 mm, and L = 4 mm. Moreover, S_b is set to 1 for the calculation of ε̅_xx and -1 for μ̅_yy. Apparently, both ε̅_xx and μ̅_yy demonstrate a resonant behavior at 12.33 GHz and 14.76 GHz, respectively, which lead to a sufficiently wide frequency band of concurrently negative material parameters.

 

Fig. 4. Analytically derived dispersion diagram for a planar DNG slab with a width of d = 15 mm, composed of broadside SRRs on constant y-planes. The geometrical parameters of the SRRs are: w = 3 mm, s = 0.25 mm, g = 0.5 mm, h = 0.5 mm, p_x = p_y = 5 mm, and L = 5 mm. Magnetic resonance f _m0 and magnetic plasma f_mp frequency specify the range of the resulting band, while shaded regions correspond to the different k′_z cases. Blue lines describe the variation of symmetrical modes and red lines that of antisymmetrical ones.

Download Full Size | PDF

The proposed technique is, firstly, applied to a planar slab with a width of d = 15 mm, consisting of periodically-located broadside SRRs on constant y-planes. Preserving the dimensions of the configuration examined in Fig. 3(a), the complete dispersion diagram for the numerically computed parameters of the particular SRR scenario (see inlet sketch) is depicted in Fig. 4. In accordance with the analysis of the preceding paragraphs, the structure exhibits a resonant μ̅_yy and a non-resonant, smoothly varying, ε̅_xx parameter. So, by letting μ̅_zz = 1 in the expression of k′_z,, one receives k′_z = {μ̅_yy[(ω/c)² ε̅_xx - k_y ²]}^1/2, which can offer useful insight to the nature of the supported modes. It is straightforward to discern that k′_z is written as the product of the (μ̅_yy)^1/2 term, involving just the magnetic activity of the slab, and the [(ω)/c)² ε̅_xx - k_y ²]^1/2 term, representing the wavenumber’s z-component for propagation in a dielectric medium with an electric permittivity of ε̅_xx. As a consequence, it is rational to anticipate that surface modes could be separated to those which arise due to the negative permeability and those existing as an upshot of the non-unity permittivity, like the surface waves guided by a dielectric slab. For our study to be comprehensive, let us concentrate on the shaded regions of Fig. 4 and discern the subsequent cases for k′_z. Expressly for μ̅_yy < 0, if k_y < (ω/c)(ε̅_xx)^1/2, k′_z is imaginary (region I) and the field becomes evanescent (surface polariton) inside the slab, whereas (1) has only one solution. Notice that magnetic resonance f _m0 and magnetic plasma f_mp frequency determine the lower and upper limit of the band, so created. On the other hand, when k_y > (ω/c)(ε̅_xx)^1/2, k′_z is real (region II) and the field in the metamaterial turns into a standing wave, so supporting an infinite number of solutions for (1). Finally, if μ̅_yy > 0, the arrangement behaves like a dielectric slab which excites modes in the area between the light cone k_y = ω/c and the curve k_y = (ω/c)(ε̅_xx)^1/2 (region III). Analogous numerical outcomes and interpretations, with a suitable frequency scaling, are derived for a structure comprising edge-coupled SRRs, as well.

 

Fig. 5. Analytically derived dispersion diagram for a planar DNG slab with a width of d = 15 mm, composed of S-rings on constant y-planes. The geometrical parameters of the S-rings are: w = 0.5 mm, s = 0.25 mm, h = 0.5 mm, ℓ₁ = ℓ₂ = 3 mm, p_x = 5.5 mm, p_y = 2.5 mm, and L = 4 mm. Electric f _e0 and magnetic f _m0 resonance frequencies indicate the range of the resulting bandgap. Blue lines describe the variation of symmetrical modes and red lines that of antisymmetrical ones.

Download Full Size | PDF

Next, we proceed to the investigation of the surface waves guided along a planar metamaterial slab, of d = 15 mm, formed by thin wires and SRRs or S-rings. Herein, k′_z is given by the same expression as in the previous application. Nevertheless, ε̅_xx is not greater than zero any more, having a zero resonant frequency (Drude model) for the thin-wire/SRR medium and a nonzero resonant frequency (Lorentz model) for the S-ring substance. Therefore, for μ̅_yy < 0 and ε̅_xx < 0, k′_z is real and (1) gives rise to an infinite number of solutions. On the contrary, when μ̅_yy > 0 and ε̅_xx < 0 (in such structures f_mp is always smaller than the electric plasma frequency f_ep, while f _m0 is larger than the electric resonance frequency f _e0), k′_z becomes imaginary and no solution exists for (1). Remember that in order to solve (1) for imaginary k′_z values, μ̅_yy should be negative. Furthermore, as far as the S-ring configuration is concerned, for frequencies below the electric resonance, both constitutive parameters are positive and the slab supports surface modes similar to those encountered in region III of Fig. 4. All these observations are summarized in the dispersion diagram of Fig. 5, which corresponds to the numerically extracted parameters of the S-ring setup with the same dimensions as in the case of Fig. 3(b). A notable feature of the diagram is the bandgap created between resonant frequencies f _e0 and f _m0.

A noteworthy variant to the arrangements analyzed above ensues if additional SRRs or S-rings are placed on constant z-planes, so as to achieve a resonant magnetic response for the z-directed magnetic field components. Then, the μ̅_yy = μ̅_zz convention can be employed with a very satisfactory approximation, to simplify the k′_z expression to the more expedient one k′_z = [(ω/c)² μ̅_yy ε̅_xx - k_y ²]^1/2. Since the surface modes in such planar metamaterial slabs have been already thoroughly examined by several authors [24]-[26], no further details will be provided in the present study. Just for illustration, consider that when ε̅_xx > 0, only a symmetrical and an antisymmetrical surface polariton are excited in the negative permeability region. Moreover, albeit ε̅_xx increases with frequency, it is possible to obtain surface modes at higher frequencies where μ̅_yy > 0, which are of the same nature as those evolving at conventional dielectric slabs.

Overall, it is stressed that our analysis, presented till now, assumes that all numerically computed local average parameters stand for the homogeneous local constitutive parameters of the slab’s DNG medium. While this approach yields very accurate results for small k_y, it becomes incorrect as k_y increases, owing to the nonlocal effects which dominate the properties of the structure. However, it is still fairly feasible to draw useful interpretations and facilitate the comprehension of the underlying physical mechanisms that influence the dispersion characteristics of surface polaritons. To certify these issues, the subsequent section conducts an in-depth investigation of periodic metamaterial setups by means of enhanced FDTD simulations, and compares the outcomes with those of the technique hitherto presented.

3. Numerical verification and results

An efficient way to evaluate the frequencies of the supported surface modes at SRR or S-ring based DNG slabs is to model the unit cell and the corresponding periodic boundary conditions in the context of the FDTD method and excite the computational domain with a broadband temporal pulse. Then, the desired frequencies are extracted from the peaks of the electromagnetic field’s Fourier spectrum and are deemed to be exact, since no specific prerequisite about the form of the field has been imposed. Conversely, numerical errors due to the discretization of the continuous space are easily controlled by, suitably, reducing the spatial step of the grid. To further support this statement, Table 1 presents the mode frequencies of the S-ring arrangement with the geometrical parameters of Fig. 5, for two different spatial discretization steps, namely λ _min/60 and λ _min/120, with λ _min the wavelength at the maximum frequency of 20 GHz. Simulations are performed for q_y = 0.45π/p_y, where the Bloch number q_y plays the role of k_y in analytical calculations, simply indicating its use as the y-directed wavenumber component. Nevertheless, contrary to k_y, q_y must belong to the first Brillouin zone, that is −π/p_y < q_y < π/p_y. So, the specific q_y value employed for the derivation of Table 1 lies in the middle of the prior zone, rendering the particular comparisons representative of all numerical simulations presented hereafter. The differences between the frequencies are in the order of 1 decimal point, which is considered very satisfactory for the present analysis. Concerning rounding errors, they do not have any effect on the results, since in every simulation double precision arithmetic (14 decimal digits) has been utilized, while field values are of the 10⁰ order. Furthermore, the only converging procedure involved in the method is related to the time needed for the steady state to establish. By comparing the Fourier spectrum for different simulation times it has been found that 10⁴ steps are enough for the transient phenomena to decay and accurate results to be obtained. Having verified the accuracy of the numerical scheme, it can be, now, implemented to derive the actual dispersion diagrams of metamaterial arrangements.

Table 1. Mode frequencies for the S-ring arrangement of Fig. 5 and q_y = 0.45π/p_y.

View Table

 

Fig. 6. Analytically (solid lines) and numerically (dots) derived dispersion curves for a planar DNG slab with a width of d = 15 mm, composed of broadside SRRs on constant y-planes. The geometrical parameters of the SRRs are: w = 3 mm, s = 0.25 mm, g = 0.5 mm, h = 0.5 mm, p_x = p_y = 5 mm, and L = 5 mm.

Download Full Size | PDF

Let us start our analysis with a slab composed of broadside SRRs placed on constant z-planes, which, for the sake of comparison, has the same dimensions as the one in Fig. 4. Figure 6 displays the dispersion curves for a set of the supported modes, calculated both numerically (red dots) and via the analytical methodology of the previous section (solid lines). Bearing in mind that the unit cell is, almost, symmetric on both sides of the y-plane which divides it into two equal parts, the dispersion curves for -π/p_y < q_y < 0 are expected to be axially symmetric with those for 0 < q_y < π/p_y. Thus, only q_y values between 0 and π/p_y are considered hereafter. Obviously, all modes predicted through the theoretical analysis are, also, computed by the FDTD simulations. Particularly, there is a very good quantitative agreement between the analytical and numerical curves for frequencies in the DNG band and below it. Conversely, in the region beyond 12.37 GHz and for q_y > 0.6π/p_y, that is for q_y approaching the upper bound of the Brillouin zone, the discrepancies between the two processes become more intense. The explanation for this behavior rests on the fact that the analytical curves of the previous modes grow as q_y increases, in contrast to the respective numerical ones, which exhibit a zero inclination for q_y near π/p_y. Notice that, by definition, the dispersion curves of every periodic structure have a zero inclination at the end of the Brillouin zone. On the other hand, the analytical curves remain practically stable in the range of 11.68-12.37 GHz, like the numerical ones, justifying the satisfactory agreement between the corresponding results in the previous range.

 

Fig. 7. Numerically derived dispersion diagram for a planar DNG slab with a width of d = 15 mm, composed of S-rings on constant y-planes. The geometrical parameters of the S-rings are: w = 0.5 mm, s = 0.25 mm, h = 0.5 mm, ℓ₁ = ℓ₂ = 3 mm, p_x = 5.5 mm, p_y = 2.5 mm, and L = 4 mm.

Download Full Size | PDF

Moving to the second application, Fig. 7 shows the dispersion diagram of an S-ring based slab with the geometrical parameters of the setup explored in Fig. 5. From the results, one can clearly distinguish the bandgap between the numerically evaluated electric and magnetic resonance frequencies f′_e0 and f _m0, respectively, predicted, also, by our theoretical formulation. Similarly, the shape and spectral distribution of surface plasmons excited above f′_m0 and below f′_e0 - as well as everywhere in the entire frequency range - are successfully extracted. Again, it is interesting to observe that all numerical dispersion curves attain a zero inclination angle for q_y = π/p_y. Further comparison of the dispersion diagrams in Figs 5 and 7, indicates some nontrivial discrepancies in the values of resonant frequencies for constant q_y, especially in the band above f′_m0. As a matter of fact, frequencies f′_e0, f′_m0 are shifted to somewhat larger values with respect to the theoretically extracted f _e0, f _m0 ones. However, such a behavior does not, by any means, degrade the expediency and versatility of the analytical approach as a qualitative and computational tool for planar metamaterial arrangements.

The final example studies a thin slab with a width of d = 5 mm, that consists of one broadside SRR toward the z direction. Basically, this is a case of crucial importance, since one of the most remarkable properties of SNG slabs is their competence to, effectively, guide the electromagnetic energy, no matter how small their width is. Therefore, Fig. 8 presents its dispersion diagram maintaining the dimensions of the SRR used in Fig. 6. It is obvious that surface waves exist in all three regions of the diagram, described in the previous section. Actually, the numerically derived dispersion curves are in very good agreement, concerning their shape, with the analytical ones, thus proving the accuracy of the homogenized model based method. Furthermore, as in the case of Fig. 6, for frequencies above the magnetic plasma one, there is a good quantitative agreement when Bloch number q_y < 0.6π/p_y. Compared to the outcomes of Fig. 6, the most prominent difference is that for small slab widths, the number of the excited modes decreases. This worth mentioning phenomenon is encountered in the analytical model too, where a sparser modal distribution in the frequency range is accomplished.

 

Fig. 8. Analytically (solid lines) and numerically (dots) derived dispersion curves for a thin DNG slab with a width of d = 5 mm, composed of one broadside SRR in the z-direction. The parameters of the SRR are: w = 3 mm, s = 0.25 mm, g = 0.5 mm, and h = 0.5 mm.

Download Full Size | PDF

4. Conclusion

The dispersion characteristics of the surface waves which propagate parallel to the edges of planar slabs composed of realistic DNG or SNG media, have been systematically examined in this paper through an efficient analytical and numerical methodology. Based on the solution of the theoretical dispersion relation, under the presumption that the metamaterial is local and homogeneous, the former scheme uses the average permittivity and permeability of the unit cell’s artificial substance as the analysis’ constitutive parameters. Also, to considerably advance the overall accuracy, these average parameters are computationally determined, instead of being extracted via existing analytical expressions. In contrast, the numerical formulation is founded on a precise FDTD technique combined with periodic boundary conditions. To assert these merits, several DNG or SNG structures are explored and elaborate comparisons are conducted.

Indisputably, the numerical approach provides more accurate results for any value of the transverse wavenumber component, since no extra convention about the field is made. Potential error sources stem from the spatial discretization and they can be seriously suppressed via finer grid resolutions, at the expense of increased overheads. Conversely, the analytical algorithm is easy to implement, although it can lead to some inaccuracies near the first Brillouin zone. This issue is chiefly accredited to the locality assumption for the field embodied in the formulation. To overcome this difficulty, nonlocal effects must be, somehow, included in the surface modes dispersion equation, possibly by altering its basic form. It is, finally, stated that, even when the average parameters are evaluated by a single numerical simulation, they can be safely used in the theoretical formulae to model different arrangements with the same unit cell, a fact which proves the robustness of the homogenized model technique over the numerical approach.