© M. Naeem & Y. Hao, Published by EDP Sciences, 2016

1 Introduction

Recent advances in electromagntics applications requiring non-natural material properties, such as tranformation optics [1] and field transformation [2], have accentuated the importance of materials design. Luneberg lens and its variant Maxwell’s fish-eye lens, for instance, require continuously varying material properties, namely, their refractive index varies smoothly along the radial direction from high (in the centre) to low (at the edge). One possible practical realisation of Maxwell’s fisheye is the use of metamaterials, that is, using metallic elements that mimic the required dielectric properties by varying the impedance of two parallel metallic plates. A rather straightforward avenue to vary the impedance of a parallel plate is to radially vary the separation between the plates [3, 4]. The main disadvantage of using metallic components to design the effective properties is the material loss in the conducting elements. Alternatively, a material with the desired electromagnetic properties can be designed, in principle, by varying the volume fraction of the inclusion to host materials, or by altering their dielectric contrast [5]. The effective permittivity of the resultant material remains in between the minimum and maximum values of the inclusion and host materials depending on the following bounds:(1a)

(1b)where v₁ … v_n − 1 are the volume fractions (inclusions to total volume of the composite) of ε₁ … ε_n − 1, respectively; and v_n = v_total − (v₁+ … + v_n−1) is the volume fraction of the host with permittivity ε_n. It is worth noting that these bounds have limited validity, that is, in the quasi-static regime only. In this context, only two off-the-shelf materials can be utilized to span a range of dielectric or magnetic properties. As a generalization, the inclusion and host materials can be considered as two elementary building blocks, known as metamaterial bits, and a spatial arrangement of these bits, dubbed a metamaterial byte, results in homegenized properties that can be judiciously manipulated to any permittivity value in between the Wiener bounds [6]. The design procedure, however, requires improved material modelling schemes that can be efficiently used for material homogenization. The current homogenization schemes have their own pros and cons; predominately, the homogenization schemes are formulated based on the analytical treatments for given arrangement of circumstances only. For instance, an infinitely large periodic structure can be homogenized, below the diffraction limit, by using the powerful tools of the Bloch analysis originally developed in solid state physics. Though this procedure provides accurate and numerically efficient solution to infinite many elementary particles in a crystal, it inherently ignores the edge effect, which is often desired in solving mesoscale electromagnetic problems. Here, we introduce a point-dipole based model [7] for the homogenization of the periodic or aperiodic composites. In its formulation, the volumetric polarization currents inside the dielectric body are replaced by the equivalent dipole moments whose weights are then determined by enforcing a consistency condition for the electric fields inside the dielectric body. As a result, the proposed method of equivalent dipole moments (MEDM) does not suffer from the singularities encountered in the conventional method of moment formulations that are based on the Green’s function approaches [8], where a superposition is made of the fields generated by infinitely small concentrated current sources.

2 Limitations of existing point-dipole homogenization models

2.1 Clausius-Mossoti formulation

Consider an infinitely large periodic arrangement of scatterers (having equal axial-ratio) with relative permittivity ε_r embedded in a background material of permittivity ε₀ as depicted in Figure 1. Clausius-Mossoti (CM) derivation assumes a fictitious cavity at the interior of a heterogeneous material: the contribution of the scatterers exterior to the cavity is calculated by integrating the contribution of the surface charges at the cavity, while the field of the interior scatterers is added one by one. The local (loc) field at the fictitious cavity is then the sum of the incident (i) field, and the field contribution by the rest of the scatterers in the background (test scatterer removed) as:(2)

Figure 1. Periodic arrangement of dielectric scatterers of constitutive parameters (μ0,ε0) embedded in a medium of background medium of permittivity (μ0, εr). Arrows inside the scatterers indicate the depolarization field under the application of an incident field Ei. The symbol ECell represents the total field in a cubic unit-cell.

In the case of identical inclusions the field observed at the centre of the removed scatter, E_loc, is the sum of the fields scattered by the rest of the scatterers in the array; its value is equal to the field inside a standalone scatterer [9, 10]. The local field can thus be given as:(3)

The employment of the above expression for the local field in the polarization formula of the test scatterer P = n₀αE_loc yields:(4)where, n₀ is the number of scatterers per unit volume, and α is the polarizability of the inclusion particles. The above equation can be rearranged to give:(5)

By equating equation (5) with P = ε₀(ε_r − 1)E_i, we get:(6)

The above equation can further be rearranged to give a closed-form expression for the effective permittivity of the heterogeneous material; expressed in terms of the polarizabilities of the scatterers, the Clausius-Mossotti equation [10] is then obtained as:(7)

Notice that the CM equation suffers from the singularity for n₀α/3ε₀ = 1, which is called the Mossotti catastrophe, which can be dealt with by an accurate representation of the cavity field, but the treatment would be beyond the scope of this investigation.

The difference between the Maxwell-Garnett (MG) [11] and CM is that the MG writes the effective permittivity of the heterogeneous material in terms of the constituent parameters (instead of polarizabilities) as:(8)where f = n₀v₀ represents the volume fraction of inclusion to host, v₀ is the volume of an inclusion scatterer. Many extensions to the classical MG model have been made since its inception to derive a closed-form expression for the non-spherical in order and multiphase inclusions [12, 13].

2.2 Brief description of classic homogenization techniques

We briefly review other classical techniques that we will use for comparison with the proposed scheme in Section 4. One widely known permittivity extraction technique is Nicolson-Ross-Weir (NRW) method, which involves the calculation of the scattering parameters from the test sample; except for a few analytical solutions of well-defined shapes for the reflection and transmission coefficients, the extraction is usually performed using numerical calculations. Effective parameters are then calculated using the closed-form expressions given in [14]. However, the extraction of the effective parameters when the sample thickness is close to the wavelength (resulting in the Bragg resonances) has been a matter of concern while using the NRW technique. We use a commercial software (CST Studio Suite) which uses an improved homogenization model based on the reflection and transmission coefficients [15] and performs accurately even close to the Bragg resonances. The two scale convergence homogenization (TSCH) technique (discussed in detail in [16]) is based on the numerical resolution of an annex static problem to compute the effective permittivity of a periodic medium with arbitrary shaped inclusions.

3 Proposed homogenization scheme

Most of the classical homogenization techniques have one commonality, that is, the effective parameters they extract are based on infinitely large arrangement of periodic scatterers, where we either ignore or calculate the approximate interaction field between the scatterers. Practical problems, however, mostly involve finite size geometries. We propose to solve a finite size array of custom-arranged scatterers by replacing each inclusion scatterer in the composite by its equivalent dipole moment for the scattered field representation outside the scatterer, and performing the analytical treatment of the field scattered inside the unit-cell. Unlike the conventional surface-based discretization using the RWG (Rao-Wilton-Glisson) basis in the full-wave electromagnetic models, the proposed method considers each scatterer as a basis function. The mutual interaction between the scatterers and the self impedance can thus be calculated analytically. The numerical computation involved is minor, that is, inversion of the coupling matrix. Once the scattered field is known, the effective parameter extraction involves the homogenization of the field in the considered unit-cell and homogenization of the equivalent current in the test scatterer. By the application of the field consistency condition, we have(9)where is the volume of the test scatterer. Field consistency condition allows for calculating the effective values of the permittivity of the material.(10)where Pol ∈ {x, y, z} is the polarization, and the average local field is calculated in a cubic unit-cell of the edge-length which is equal to the period, P, of the lattice. For a centersymmetric scatterer which is placed at the origin of the rectangular coordinate system, we can calculate the total field in the cell as:(11a)

(11b)

(11c)

The cell-field includes in part the Rayleigh scattering field, and in part the dipole radiation field. The optimal calculations for finite size geometries thus require the full-wave modelling treatment. The equivalent current in the definition of equation (10) also includes the mutual interactions among the particles. By definition, this field is constant inside the test scatterer and zero elsewhere, which allows us to integrate the equivalent current inside a basis of radius r only.(12a)

(12b)

(12c)

Notice that the symbols and represent the field and current in a unit-cell volume, that is, in the units of volts meter² and amperes meter, respectively. Analytical considerations of the proposed homogenization procedure (i.e., Eq. (10)) reveal that the averaging implicitly caters for the phase of the scattered field as well as the equivalent current; in consequence, while the proposed method inherently sidesteps the low frequency breakdown problem, it also allows for the calculation of the effective permittivity at considerably high frequencies.

The effective permeability of the composite can be calculated by duality, but the proposed scheme need to be effectively adapted to cater for the electric and magnetic dipole coupling to include both the electric and magnetic dipole moments simultaneously. Moreover, the dynamic expressions for the self-term and the dipole scattered fields can efficiently calculate the dipole moments (DMs) of both the dielectric and metallic scatterers. We justify this assertion by calculating the dipole moments (using the proposed scheme [7]) for a standalone dielectric and metallic scatterer, the results are as follows:(13a)

(13b)

In the above equations, is the volume of the scatterer, and E_i represents the amplitude of the incident field. The calculated values are consistent with those provided in the Eqs. 6.110 and 6.115 of the reference [17]; note that the equation (13b) has been evaluated by the application of high conductivity limit.

4 Results

We considered a d-periodic arrangement of spherical scatterers in a three-dimensional free space (relative permittivity = 1) with the aim of establishing the relation or identifying the differences among various homogenization schemes (Figure 2).

Figure 2. Effective permittivity extraction of a periodic material using various homogenization techniques. We used full wave-finite element method based Nicolson Ross Weir (FEM-NRW) technique, Maxwell-Garnett (MG) periodic, MG finite along one dimension (z-direction), MG finite (finite along each direction), TSCH, and the proposed technique (MEDM) in order to extract the effective parameters. A cubic sample of the composite is considered, and the full geometric and simulation settings are given in Table 1. (a) Effective permittivity as a function of the permittivity of the inclusion material. The lattice period d is three times the edge-length of the basis a. (b) Effective permittivity as a function of the lattice period. The relative permittivity of the inclusion is 6.

We have used the Maxwell-Garnett technique for a fully periodic medium (unit-cell, volume fraction f = a ³/d ³), finite width h = (N − 1)d + a in z direction (f = Na ³/hd ²), and finite size in three dimensions (f = (Na)³/h ³). The proposed scheme involved a finite volume of periodic scatterers in three orthogonal dimensions and considered a unit-cell in the center of the array to include a considerable amount of interaction with the neighbouring scatterers. The effective permittivity calculations (Eq. (10)) thus resulted in close agreement to that of the existing schemes, except close to high volume fraction of inclusion to host (Figure 2b), where problem is the poor approximation of the local field (Eq. (3)). Bearing in mind the electrostatic nature of the MG, we compared it to the proposed homogenization scheme and found promising results even for smaller lattice constant and high dielectric contrast between the inclusion and background media. Though the TSCH scheme has the ability to discretize arbitrary shaped objects in the periodic boundary condition settings, we found it not as accurate as the MG or MEDM.

For numerical homogenization, we used periodic boundary conditions in two orthogonal axis that are perpendicular to the plane of incidence, and used finite size array (5-layers) in the plane of incidence. The results can be discerned from the MEDM due to the finite size of the array in the incidence plane, which treats the material as a film rather than a bulk. To explain this in more detail, we considered a single layer of scatterers periodic in the other two dimensions, we observed that the effective permittivity is higher in the case of sheet and approaches to that of the bulk permittivity when the interaction from the neighbouring layers is considered, as shown in Figure 3. The results are consistent to that of the published data which claims that the effective parameter values are higher in sheets and reduce to that of bulk parameter values when the contribution from the neighbouring layers is considered [18]. A commercial electromagnetic software (CST) has been used in frequency domain in order to extract the effective parameters. The CST simulation time (21 s) was comparable to that of the MEDM (19 s), which was implemented in MATLAB. We will demonstrate in a future communication that the simulation time improves dramatically in simulating inhomogeneous composites such as composite gradients, where the commercial software can not leverage the computations by applying the periodic boundary conditions.

Figure 3. The effect of the finite thickness of the sample on the extracted parameters. The extracted permittivity of two-dimensional arrangement of scatterers (in good agreement for both the CST and proposed scheme) is higher than the bulk permittivity. Full geometric and simulation settings are the same as Figure 2b except for the single layer consideration.

5 Conclusion

Due to its rigorous formulation and promising accuracy, the proposed technique allows for the numerical homogenization of finite-sized composites at a reduced computational cost as compared to the conventional full-wave methods. Numerical calculations reveal that the effective permittivity strongly depends on the volume fraction and dielectric contrast (inclusion to host) of the unit-cell considered; the lower the volume fraction or dielectric contrast, the lower is the permittivity, and vice versa. It would be fruitful to extend the research and make the scattering analysis feasible close to the resonances in the scatterers, and due to the aperiodic arrangement of scatterers.

Acknowledgments

The authors would like to thank Dr. Benjamin Vial for providing the TSCH results for comparison with the proposed method.

References

All Tables

All Figures

	Figure 1. Periodic arrangement of dielectric scatterers of constitutive parameters (μ0,ε0) embedded in a medium of background medium of permittivity (μ0, εr). Arrows inside the scatterers indicate the depolarization field under the application of an incident field Ei. The symbol ECell represents the total field in a cubic unit-cell.
In the text

Figure 2. Effective permittivity extraction of a periodic material using various homogenization techniques. We used full wave-finite element method based Nicolson Ross Weir (FEM-NRW) technique, Maxwell-Garnett (MG) periodic, MG finite along one dimension (z-direction), MG finite (finite along each direction), TSCH, and the proposed technique (MEDM) in order to extract the effective parameters. A cubic sample of the composite is considered, and the full geometric and simulation settings are given in Table 1. (a) Effective permittivity as a function of the permittivity of the inclusion material. The lattice period d is three times the edge-length of the basis a. (b) Effective permittivity as a function of the lattice period. The relative permittivity of the inclusion is 6.

In the text

	Figure 3. The effect of the finite thickness of the sample on the extracted parameters. The extracted permittivity of two-dimensional arrangement of scatterers (in good agreement for both the CST and proposed scheme) is higher than the bulk permittivity. Full geometric and simulation settings are the same as Figure 2b except for the single layer consideration.
In the text