https://doi.org/10.1051/epjam/2016004
Research Article
Homogenization of composites using fullwave pointdipole model
School of Electronic Engineering and Computer Science, Queen Mary University of London, Mile End Road, London
E1 4NS, UK
^{*} email: m.naeem@qmul.ac.uk
Received:
18
April
2016
Accepted:
30
May
2016
Published online: 26 July 2016
We propose a fullwave pointdipolebased scheme for accurate material homogenization of composites. Unlike conventional homogenization approaches which employ periodic boundary conditions for the effective permittivity extraction, the proposed technique can efficiently compute the interactions between the scatterers even in inhomogeneous mixtures as well as aperiodic arrangement of a finite number of scatterers. The method utilizes each inclusion as a basis function, and thereby sidesteps the computational burden to discretize each inclusion scatterer. In addition, the supplementary analytical calculations provide closedform expressions to deal with the field singularities and to calculate the mutual interaction between the scatterers. We demonstrate that the proposed homogenization scheme can precisely extract the composite permittivity profile of a finite size specimen.
Key words: Material homogenization / Field transformation / Transformation optics
© M. Naeem & Y. Hao, Published by EDP Sciences, 2016
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Recent advances in electromagntics applications requiring nonnatural 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 fisheye 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_{1} … v_{n − 1} are the volume fractions (inclusions to total volume of the composite) of ε_{1} … ε_{n − 1}, respectively; and v_{n} = v_{total} − (v_{1}+ … + 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 quasistatic regime only. In this context, only two offtheshelf 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 pointdipole 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.
Specifications of the simulated geometries. The superscript 3 indicates the periodicity in the three orthogonal directions. Host relative permittivity is set to one in each case.
2 Limitations of existing pointdipole homogenization models
2.1 ClausiusMossoti formulation
Consider an infinitely large periodic arrangement of scatterers (having equal axialratio) with relative permittivity ε_{r} embedded in a background material of permittivity ε_{0} as depicted in Figure 1. ClausiusMossoti (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 E_{i}. The symbol E_{Cell} represents the total field in a cubic unitcell. 
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_{0}αE_{loc} yields:(4)where, n_{0} 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 = ε_{0}(ε_{r} − 1)E_{i}, we get:(6)
The above equation can further be rearranged to give a closedform expression for the effective permittivity of the heterogeneous material; expressed in terms of the polarizabilities of the scatterers, the ClausiusMossotti equation [10] is then obtained as:(7)
Notice that the CM equation suffers from the singularity for n_{0}α/3ε_{0} = 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 MaxwellGarnett (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_{0}v_{0} represents the volume fraction of inclusion to host, v_{0} is the volume of an inclusion scatterer. Many extensions to the classical MG model have been made since its inception to derive a closedform expression for the nonspherical 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 NicolsonRossWeir (NRW) method, which involves the calculation of the scattering parameters from the test sample; except for a few analytical solutions of welldefined shapes for the reflection and transmission coefficients, the extraction is usually performed using numerical calculations. Effective parameters are then calculated using the closedform 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 customarranged 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 unitcell. Unlike the conventional surfacebased discretization using the RWG (RaoWiltonGlisson) basis in the fullwave 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 unitcell 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 unitcell of the edgelength 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)
The cellfield includes in part the Rayleigh scattering field, and in part the dipole radiation field. The optimal calculations for finite size geometries thus require the fullwave 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)
Notice that the symbols and represent the field and current in a unitcell volume, that is, in the units of volts meter^{2} 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 selfterm 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)
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 dperiodic arrangement of spherical scatterers in a threedimensional 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 wavefinite element method based Nicolson Ross Weir (FEMNRW) technique, MaxwellGarnett (MG) periodic, MG finite along one dimension (zdirection), 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 edgelength 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 MaxwellGarnett technique for a fully periodic medium (unitcell, volume fraction f = a ^{3}/d ^{3}), finite width h = (N − 1)d + a in z direction (f = Na ^{3}/hd ^{2}), and finite size in three dimensions (f = (Na)^{3}/h ^{3}). The proposed scheme involved a finite volume of periodic scatterers in three orthogonal dimensions and considered a unitcell 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 (5layers) 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 twodimensional 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 finitesized composites at a reduced computational cost as compared to the conventional fullwave methods. Numerical calculations reveal that the effective permittivity strongly depends on the volume fraction and dielectric contrast (inclusion to host) of the unitcell 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
 J.B. Pendry, D. Schurig, D.R. Smith, Controlling electromagnetic fields, Science 312 (2006) 1780–1782. [CrossRef] [MathSciNet] [PubMed] (In the text)
 S. Jain, M. Abdelmageed, R. Mittra, Flatlens design using field transformation and its comparison with those based on transformation optics and ray optics, IEEE Antennas and Wireless Propagation Letters 12 (2013) 4–7. [CrossRef] (In the text)
 J. Liu, R. Mendis, D.M. Mittleman, A Maxwell’s fish eye lens for the terahertz region, Applied Physics Letters 103 (2013) 031104. [CrossRef] (In the text)
 R. Mendis, D.M. Mittleman, Artificial dielectrics, IEEE Microwave Magazine 15 (2014) 34–42. [CrossRef] (In the text)
 S. Kubo, A. Diaz, Y. Tang, T.S. Mayer, I.C. Khoo, T.E. Mallouk, Tunability of the refractive index of gold nanoparticle dispersions, Nano Letters 7 (2007) 3418–3423. [CrossRef] (In the text)
 C. Della Giovampaola, N. Engheta, Digital metamaterials, Nature Materials 13 (2014) 1115–1121. [CrossRef] (In the text)
 M. Naeem, Scattering and Absorption Analysis of Radomes Using the Method of Equivalent Dipole Moments, PhD Thesis, Chalmers University of Technology, 2011. (In the text)
 R.D. Graglia, On the numerical integration of the linear shape functions times the 3D Green’s function or its gradient on a plane triangle, IEEE Transactions on Antennas and Propagation 41 (1993) 1448–1455. [CrossRef] (In the text)
 D.E. Aspnes, Localfield effects and effectivemedium theory: a microscopic perspective, American Journal of Physics 50 (1982) 704. [NASA ADS] [CrossRef] (In the text)
 L. Tsang, Scattering of Electromagnetic Waves, Wiley, New York, 2000. [CrossRef] (In the text)
 J.C. Maxwell Garnett, Colours in metal glasses and in metallic films, Transactions of the Royal Society 203 (1904) 385–420. [NASA ADS] [CrossRef] (In the text)
 S. Bosch, J. FerréBorrull, N. Leinfellner, A. Canillas, Effective dielectric function of mixtures of three or more materials: a numerical procedure for computations, Surface Science 453 (2000) 9–17. [NASA ADS] [CrossRef] (In the text)
 M.Y. Koledintseva, R.E. Dubroff, R.W. Schwartz, Mixtures containing conducting particles at optical frequencies, Progress In Electromagnetics Research 63 (2006) 223–242. [CrossRef] (In the text)
 A.M. Nicolson, G.F. Ross, Measurement of the intrinsic properties of materials by timedomain techniques, IEEE Transactions on Instrumentation and Measurement 19 (1970) 377–382. [CrossRef] (In the text)
 X. Chen, T.M. Grzegorczyk, B.I. Wu, J. Pacheco, J.A. Kong, Robust method to retrieve the constitutive effective parameters of metamaterials, Physical Review E 70 (2004) 016608. [CrossRef] (In the text)
 S. Guenneau, F. Zolla, Homogenization of threedimensional finite photonic crystals, Progress in Electromagnetics Research 27 (2000) 91–127. [CrossRef] (In the text)
 R.F. Harrington, TimeHarmonic Electromagnetic Fields, IEEE Press, New Jersey, 2001. [CrossRef] (In the text)
 C.A. Grimes, The effective permeability of granular thin films, IEEE Transactions on Magnetics 29 (1993) 4092–4094. [CrossRef] (In the text)
Cite this article as: Naeem M & Hao Y: Homogenization of composites using fullwave pointdipole model. EPJ Appl. Metamat. 2016, 3, 6.
All Tables
Specifications of the simulated geometries. The superscript 3 indicates the periodicity in the three orthogonal directions. Host relative permittivity is set to one in each case.
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 E_{i}. The symbol E_{Cell} represents the total field in a cubic unitcell. 

In the text 
Figure 2.
Effective permittivity extraction of a periodic material using various homogenization techniques. We used full wavefinite element method based Nicolson Ross Weir (FEMNRW) technique, MaxwellGarnett (MG) periodic, MG finite along one dimension (zdirection), 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 edgelength 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 twodimensional 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 