Computational electromagnetics and the rational design of new dielectric heterostructures

Progress in Materials Science 48 (2003) 373–456 www.elsevier.com/locate/pmatsci

Computational electromagnetics and the rational design of new dielectric heterostructures C. Brosseaua,b,
, A. Beroualc Laboratoire d’E´lectronique et Syste`mes de Te´le´communications, Universite´ de Bretagne Occidentale, B.P. 809, 6 avenue Le Gorgeu, 29285 Brest Cedex, France b ´partement De de Physique, Universite´ de Bretagne Occidentale, Brest, France c Centre de Ge´nie E´lectrique de Lyon, E´cole Centrale de Lyon, B.P. 163, 36 avenue Guy de Collongue, 69131 Ecully Cedex, France a

Accepted 1 July 2002

Abstract Dielectric properties of heterogeneous materials for various condensed-matter systems have been gaining world-wide attention over the past 50 or so years in the design (or engineering) of materials structures for desired properties and functional purposes. These applications range from cable and current limiters to sensors. These multiscale systems lead to challenging problems of connecting micro- or meso-structural features to macroscopic materials response, i.e. permittivity, conductivity. This article first reviews progress made at that time of the underlying physics of dielectric heterostructures and points out the missing elements that have led to a resurgence of interest in these and related materials. Recent advances in computational electromagnetics provide unparalleled control over morphology in this class of materials to produce a seemingly unlimited number of exquisitely structured materials endowed with tailored electromagnetic, and other physical properties. In the text to follow, we illustrate how an ab initio computational technique can be used to accurately characterize structure–dielectric property relationships of periodic heterostructures in the quasistatic limit. More specifically, we have carried out two-dimensional (2D) and three-dimensional (3D) numerical studies of two-component materials in which equal-sized inclusions, with shape and orientation and possibly fused together, are fixed in a periodic square (2D) or cubic (3D) array. Boundary-integral equations (BIE) are derived from Green’s theorem and are solved for the local field with appropriate periodicity conditions on a unit cell of the structures using the field calculation package PHI3D. A number of illustrative examples shows how this computational technique can provide very accurate predictions for the complex effective permittivity




Corresponding author. Tel.: +33-2-98-01-61-05; fax +33-2-98-01-61-31. E-mail addresses: [email protected] (C. Brosseau), [email protected] (A. Beroual).

0079-6425/03/$ - see front matter # 2003 Elsevier Science Ltd. All rights reserved. PII: S0079-6425(02)00013-0

374

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

of translationally-invariant heterostructures. The performance of the method is also compared with those of other computational and analytical techniques. We comment on how this computational method helps identify some important characteristics for rationalizing and predicting the structure of composite materials in terms of the nature, size, shape and orientation of their constituents. # 2003 Elsevier Science Ltd. All rights reserved.

Contents 1. Introduction: scope of the review.................................................................................... 374 2. Background: modeling the permittivity of composites....................................................379 3. Historical perspective ......................................................................................................382 4. Numerical procedure.......................................................................................................390 4.1. Problem statement..................................................................................................391 4.2. Algorithm ...............................................................................................................393 5. Numerical examples and discussion ................................................................................ 396 5.1. Two-dimensional simulations ................................................................................. 397 5.1.1. Square lattice and 2-pole approximation ...................................................397 5.1.2. ‘‘Random’’ configurations .......................................................................... 399 5.2. Three-dimensional simulations ............................................................................... 402 5.2.1. Nonabsorbing configurations ..................................................................... 402 5.2.2. Absorbing configurations, scaling and critical exponents ..........................415 5.2.3. Multilayered structures............................................................................... 426 5.2.4. Comparison with experiment ..................................................................... 429 5.3. Discussion and final remarks.................................................................................. 435 6. Conclusions and future directions................................................................................... 438 Acknowledgements...............................................................................................................444 Appendix A ..........................................................................................................................445 Appendix B ..........................................................................................................................445 References ............................................................................................................................448

1. Introduction: scope of the review The dielectric properties of disordered heterogeneous materials [1–4] have received considerable fundamental attention since nearly the time of Maxwell and these investigations continue to be centrally important for developing a deeper understanding

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

375

leading to a wide variety of innovative technological applications in many branches of engineering. Examples include, but are far from limited to, electrical cable sheathing, electric paints for electromagnetic shielding, xerographic inks, static electricity dissipation devices, current-limiting thermistors, and sensors. Although, experiments and simulations designed to probe the dielectric behavior of composite materials containing a mixture of conducting and insulating components have uncovered quite different pieces of information, the a priori prediction of dielectric properties of these media remains an outstanding problem in materials physics. Fundamental knowledge of the polarization mechanisms in this important form of condensed matter is far from complete. Recent investigations have brought into focus that many pertinent questions remain to be answered. For example, the subtleties of the proper treatment of anisotropy and spatial inhomogeneity is still not comprehensively understood and simulation of random heterostructures is a field still in its infancy. An understanding of the relationship between the structure of the system and the macroscopic electromagnetic properties would provide insight into the design of composite materials with properties that depend on the spatial arrangement of the constituents in the mixture. Rationalizing and predicting the electromagnetic properties in terms of the nature and the structure of their constituent has, however, proved notoriously difficult for several reasons. First, numerical solution of the partial differential equation for the local electric field which occurs in the quasistatic description of random heterostructures is very cumbersome (even when the detailed features of the microstructure is known) in practice because it requires the knowledge of the pair and higher correlation functions that specify the average microscopic arrangement of the constituents. Since it is impossible to deal with all analytical specifications of these correlation functions, the statistical approach of randomness and connectedness of a real composite system is reduced on the investigation of of simple models of heterostructures [5,6]. Second, while calculations on composite materials were once the domain of empirical methods, they are now dominated by ab initio numerical aproaches. These impressive strides have been made possible by a confluence of substantial gains in hardware performance, efforts directed toward improvements in the computational economy of software, and the developments of new mathematical formalisms. The everincreasing need for devices capable of operating at high power levels, high temperatures and chemically hostile environments and also being inexpensive, lightweight, easily shaped and bonded to or embedded in a variety of surfaces, is driving researchers to develop composite materials which are capable of handling these technological challenges. In the most general sense, a composite material can be defined as the combination of many materials (or the same material in different states) present as separated components with different geometric and physical properties, and mixed together to form structures that may improve a specific property [7,8]. Industrial interest centers on creating composites with the combined properties of the components or less expensive composites by combining costly filler particles and inexpensive polymers with limited property reduction; academic interest focuses on the causes of adsorption, aggregation and so on, as well as mixing strategies. A canonical example is

376

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

carbon black/epoxy resin composite. It is by now well understood that the presence of a small concentration of conductive filler such as carbon black particles in the thermosetting polymer host can affect dramatically the dielectric properties of the system. The prototype of these remarkable phenomena is the percolation transition. At low filler concentrations, the samples reflect the properties of the matrix material. At a concentration threshold, there is a crossover between an insulating behavior (in which the carbon black particles behave as isolated aggregates) and a conductor behavior corresponding to the formation of infinitely large connected clusters of carbon black aggregates, creating a connection between the electrodes, i.e. the socalled percolation paths or channels. The formation of a conducting network depends on parameters such as the rheology of polymer melt during the mixing process, the wettability of the particles by the polymer, the details of the solidification process after mixing, and the many factors that control the mixing process, for example, time, temperature. An additional complication to the development of a detailed dielectric understanding for this system comes from the fact that the macroscopic properties depend not only on coverage, but also on how the filler surface is prepared. A thorough summary of the electromagnetic properties of filled polymer systems and other similar materials is available for the interested reader [2,9]. A wide variety of processes lead to intricate condensed-matter structures having complex morphologies with microscopic order but being macroscopically disordered. However many of them show striking similarities, for example, scaling behavior, long-range correlations. A fundamental question is whether these diverse systems are described by the same physics. In recent years some progress has been made toward answering this question by means of experiments, analytical theories and computer simulations. It goes without saying that there are a number of systems of practical interest where the answer to this question could find applications. Typical examples include filled polymers, fluid-saturated porous materials, and colloidal suspensions, to mention but a few. Each individual system is of course unique unto itself but, consistent with this diversity, there are overall similarities which one would like to explain. Numerical simulation indicates that one can in fact explain at least partially many of the observed features by simply considering a quasistatic approach, ignoring altogether the details of the microstructure, and any other complications, for example, aggregation, adsorption, and so on. This suggests that the observed dielectric behavior should be described almost completely by universal or quasiuniversal properties that are usually ascribed to structural disorder or dynamical cooperative behavior. Such behavior was given notoriety by Jonscher, who termed it ‘‘universal dielectric response’’ [10,11]. On the theoretical side there is a need for models that incorporate basic properties of individual constituents and predict the actual properties of the composites before any experimental investigations need to be performed. Newly proposed theoretical models, which are advanced with the implicit understanding that they possess general predictive powers across broad classes of heterostructures, can legitimely be scrutinized with respect to their ability to reproduce a large body of reliable experimental results as possible. Due to the increasing cost of trial-and-error approaches to determine the optimal processing parameters, industry must rely more on

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

377

numerical simulations of device procesing to find the optimal parameters. The mixing of the different components at various length scales induce changes in the microstructure which are also mirrored by a rich array of mechanical properties. The mechanical and dielectric properties of these composite materials are controlled, in large part, by the meso- and micro-structure and the interfaces between the different components, for example, the elasticity network [12]. They are materials where the primary physical properties are a consequence of the process of their assembly. For example, returning to our previous example of polymer carbon black mixtures, mechanical mixing of molten polymer chains and carbon black powder yields a morphology characterized by a supramolecular hierarchical multiscale network of polymer chains in space and time. The interatomic distances are on the microscopic level. The domain sizes and some interface widths are on the mesoscopic scale; long-range crystalline order can exist in one or more directions. The physical properties of such materials will clearly be influenced by these multiple length scales. The structure analysis in such materials is complicated due to the presence of the multitude of lengths, and also because the experimental probes used have different, and often limited, space and time resolution. This precludes reliance on simplified models invoking independent constituents of the heterostructures. Therefore, control of the sample mesostructure, and hence, macroscopic properties must be predicated upon the judicious manipulation of mixing conditions. Moreover, when dispersing a particulate conductive filler such as carbon black it is important to keep the filler volume fraction as low as posible in order to maintain the fracture toughness and the tensile properties of the polymer matrix. A recent focus of this field of research has been to correlate AFM (atomic force microscopy) observations on ultrathin sections with other physical properties, for example, local resistance, to obtain increasingly detailed information about the conduction (percolation) threshold which is evidenced in polymeric systems filled by conductive particles [4,9]. The measurements are quite difficult rendered more so by the complexity of the materials. However since the earlier experiments, researchers have continually refined them in the process answering some of the concerns about confounding factors. Over the past decade there has been plenty of progress on numerical approaches for calculating a variety of properties of heterostructures. We are developing a systematic approach to understanding these interesting physical curiosities. One component of this effort is to limit the number of ‘‘degrees of freedom’’ available to the constituents by using periodic systems, i.e. characterized by translation symmetry. To this end, the authors [13–20] and others [21–25] have described calculation techniques that use periodic models in which identical structures of constituent 1 are embedded in a crystalline fashion (Bravais lattice) in a host of constituent 2. Here, we limit our discussion to systems in the quasistatic approximation (QA), or longwavelength limit (l
). On a heuristic level, the QA can be expressed as follows: In the limit l
, where l is the wavelength of the electromagnetic radiation probing the system and
the typical scale of inhomogeneities, i.e. the correlation length for the fluctuation in the medium, a great conceptual simplification occurs because the wave cannot resolve the individual scattering centers. Consequently the medium can be said to be in a uniform average field and is characterized by an effective (average)

378

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

permittivity. Although the long-wavelength waves cannot ‘‘see’’ the inhomogeneities, the effective permittivity strongly depends on the microgeometry of the composite, i.e. the shape and topological arrangement of the components [26]. For shorter wavelengths, models of wave transport should take into account the multiple scattering losses [27–30] due to random walk diffusion of the wave from scatterer to scatterer through the material medium. In recent terms, there has been an explosion of interest in such modeling due to the development of solution techniques for partial differential equations, such as Maxwell’s equations, as well as the availability large-scale computational resources [31,32]. This article is prompted by recent progress in these endeavors. The primary purpose of this paper is to discuss an ab initio computational technique to develop interpretative and predictive descriptions of the dielectric response of periodic heterostructures in the QA. D-dimensional media will be considered (D=2 and 3 being the cases of greatest physical interest). The secondary closely connected purpose is to attempt to rationalize the relevant variables controlling the dielectric behavior of periodic heterostructures in terms of nature, size, shape, and orientation of their constituents. The informed reader will notice that the quoted literature is only a small fraction of what is published in the field. This is because the number of articles and books on the subject is by now huge, and a complete review of the significant work is beyond the purpose of this work. The layout of the paper is as follows. In the following section we attempt to outline, from the theoretical perspective, the basic physical framework in terms of which the problem of evaluating the dielectric properties of matrix-particle composites can be approached. Such composites consist of separated inclusions of material of type one embedded in a matrix (continuous phase) of material of type two. We note at the outset that this kind of composite (also termed ‘‘porphyritic’’ dispersions) is prevalent in nature. We focus on the present understanding of what we think are the most remarkable physical properties of these finely dispersed materials, observed by the state of the art experiments and interpreted with modern theories. Next for the purpose of establishing a clear context for our study, we continue by providing a brief account of the relevant historical background of some ideas and modeling efforts that are basic to our current understanding of the dielectric properties of composite materials and then highlight the role of computational electromagnetics approaches with which to evaluate the permittivity of heterostructures. With this background in mind, the implementation of the numerical approach used in this work is then set forth, which should be helpful for other physicists trying to reproduce or extend our computational results for several systems. Next, in the main part of the paper, we shall work through a number of reference two-component structures to get a feeling for scientists actively engaged in research on dielectric properties of heterostructures and also to compare our approach with others within the QA framework. Here, our computational scheme is used to investigate geometrical effects on the dielectric behavior as probed by the values of the permittivity. A direct comparison of the numerical results with our experimental results will be also discussed. Finally, the paper ends with a summary of our main results and a discussion of several open subjects that remain ahead in computational electro-

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

379

magnetics of composite materials. Two appendixes contain some explicit relations which are useful for application of the method.

2. Background: modeling the permittivity of composites Let us begin our discussion with some heuristic remarks of the basic physics of our problem. The reader may not be totally familiar with the concepts involved in modeling the dielectric behavior of heterostructures, and some simple explanations are in order. To begin with generalities, the field of dielectric properties of composite materials has come to distinguish between two different challenges. These approaches can classified as either direct or inverse. (a) Broadly speaking, in the direct problems, one tries to answer to the fundamental question: How does one calculate the dielectric macroscopic characteristics, for example, permittivity, when the properties and the arrangement of each component of the composite material are exactly known? (b) In the inverse problems, the basic question is now: How may the structure characteristics, i.e. the reconstruction of realizations of random heterogeneous media, be determined when the effective dielectric characteristics of the material are known? The first (longest-established) type of (direct) problem has led to a considerable body of phenomenology and theoretical work of varying sophistication over the years on the part of a large number of research groups worldwide, and has resulted in more or less rigorous solutions being developed for simple geometries that include spheres, cylinders and ellipsoids organized into a crystalline fashion within a matrix host. The current study falls into this class of problems. Provided that the inclusion size and mean inclusion spacing of the particulate composite material are much larger than the molecular dimensions but much smaller than the wavelength of the electromagnetic wave, the effective properties of the composite material can be calculated via homogenization theory [1,2,33,34]. There is a sizable literature on effective medium theories (EMT) to describe the macroscopic properties of random heterostructures. An EMT embodies a technique meant to bridge the gap between a detailed description of the fine-grained features of the heterostructure, and a macroscopic description which treats the medium as a completely homogeneous entity. The basic philosophy underlying most of these developments is to focus on one particular scatterer and to replace the surrounding random medium by an effective homogeneous medium. Readers familiar with the literature on statistical physics will immediately notice that homogeneization theory has much to do with mean field theory [35]. The effective medium is determined self-consistently by taking into account the fact that any other scatterer could have been chosen [1,2]. In other words, an analysis of this sort treats the homogeneity of the random medium on ‘‘large-scale’’ average. Notable theoretical papers are those of McPhedran and col-

380

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

laborators [21], Bergman [3,35–38], and Stroud [39,40], and references therein for general reviews of the field. More recent citations can be found in Ref. [2]. In spite of the increasing use of EMT over the last decades for the analysis of composite materials, progress in the assessments of EMT have been slow. Progress in research on the action of an electromagnetic wave on heterogeneous matter relies on the development of numerical techniques to solve boundary value problems. In recent times, several different strategies have been put forward for simulating the complex effective permittivity of dielectric heterostructures. First, there are finite-difference-time-domain (FDTD) methods which involve approximating the partial derivatives in the Maxwell’s equations by finite differences. In a pioneering study, Yeh [41] described the first space grid-based time integration. Taflove and Brodwin [42] demonstrated the first three-dimensional (3D) grid-based computational model of electromagnetic absorption in inhomogeneous materials. Sihvola and his group [27–29] reported the results of electromagnetic fields in dielectric heterostructures by FDTD simulations of these materials in a TEM (transverse electromagnetic) waveguide. References [31] and [43] contain detailed reviews on the FDTD methods. Spectral-domain methods have been also implemented [44,45] to obtain space and time-resolved distributions of the electromagnetic fields in inhomogeneous media. An alternative approach to these techniques is the direct resolution of boundary-integral equations (BIE), using the PHI3D field calculation package (developed at the Centre de Ge´nie Electrique de Lyon, France), to get the effective permittivity. It is now apparent that, because of its simplicity and low timing cost, this numerical approach is a computational electromagnetics methodology that is making important contributions. Prior reported work by Ghosh and Azimi [23] used a conceptually similar approach than ours. In a set of papers [13–20], we have recently described a general computational method for any geometry using a finite element algorithm and the resolution of boundaryintegrals. As we shall discuss later, the treatment can be applied to any shapes, and arbitrary distributions, however it is still limited to quasistatic situations, at the present time. Also, a major focus of our work has involved periodic systems; various parametrizations of ‘‘random’’ structures, for example, particle clustering or other local ordering, are also possible. For periodic systems, these calculations are exact. Exact means that the correct properties are computed for the given microstructure and choice of individual component properties because all internal electric multipole interactions contributing to the polarization of the material medium are taken into account [46,47]. There are strong reasons for first studying the ‘‘simplest’’ composites exhibiting the property of spatial periodicity—frequently the very essence of a property (here, the dielectric behavior) best revals itself in the simplest of cases where it occurs, unobscured by the complexity arising in more sophisticated situations. On the other hand, an exact calculation of all mutipole interactions for ‘‘random systems’’ is intractable because the complete microscopic details of the disordered configurations are not known. Numerous approximations for calculating the effective permittivity of randomly inhomogeneous materials are available in the literature [2,8,26,33,48]. The reason for this interest is the enormous variety of physical systems in which random heterogeneities occur. The means to incorporate

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

381

randomness in our lattice model has been discussed in detail previously [15]. The randomness makes the job of creation of a computer program to model dielectric properties of heterostructures a very complicated task. Computational electromagnetics of materials poses a grand challenge in systems that exhibit complex structure but can yield enormous savings in cost as well as substantially improved products. We are entering an era in which high-performance computing is coming into its own allowing true predictive simulations of complex materials to be made from information on their individual constituents and their spatial arrangement in the mixture. Despite these rapid developmental successes and quite a variety of treatments, many general problems remain that hinder further progress. Among them is the incomplete knowledge of the physical mechanisms underlying the polarization phenomena and of the disorder on the electromagnetic behavior of heterostructures. In contrast to the detailed studies of the direct problem, rather less attention and much less progress have been realized on the inverse problem. Since this case is of prime importance to experimentalists, they have resorted to applying heuristic techniques for parametrizing the dielectric response of matter [49]. A long standing goal has been the extraction of 3D structural information (size and shape of constituent particles and geometry of the associated particle boundaries) from limited experimental data, for example, using correlation functions [6,48]. Perhaps the most detailed and convincing work done to interpret the many-body nature of this problem was that of Brown [50] who showed that the effective electrical conductivity tensor of a random D-dimensional macroscopically isotropic two-component composite medium depends upon an infinite set of spatial correlation functions that statistically characterize the composite, i.e. the n-point correlation functions. Observe that correlation functions can be calculated for composites consisting of periodic arrays and random arrangements of spheres and cylinders in a matrix. An existing method of analysis, based on a truncated hierarchy of correlation functions has been also developed by several investigators, e.g. Weissberg and Prager derived expressions for the three-point correlation functions of spheres embedded randomly in a matrix [51]. Consequently, homogeneization is an issue of statistical physics [52,53]. The study of the inverse problem has not been without controversy among both theorists and experimentalists (see, for example, [54]). The success of the solution depends critically on the quality and practicality of the models. Reconstruction of the morphology of random heterogeneous media based on a record of measurements of microstructural details is a nontrivial problem with no general solution, but system-specific algorithms have been developed and demonstrated in a limited number of cases [55]. Inversion of any type of data is a notoriously difficult task. All of the inverse approaches suffer from the disadvantage of being ill-posed and often instable; the numerical instability may be overcome by regularization techniques which are well known to specialists in the field, and are discussed in detail in many places [56–59]. A key point feeding the controversy is the lack of a direct experimental method to probe in a clear cut fashion the interfaces between regions which differ in their physical characteristics [60–63]. Microphotography of thin sections of heterostructures followed by image processing can yield information about the

382

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

fractional volume of voids (porosity) and correlation functions, the statistical adequacy of which depends merely on the heterogeneity and anisotropy of the material. These data though in two-dimensional (2D), assuming invariance along the third one, can be used to adjust the reconstructed images to samples of the real heterostructures under investigation. But in practice, digitized 3D microstructures of composites are difficult to obtain and expensive. Moreover, a representative reconstruction of a medium in 3D should have the same correlation properties as those measured and averaged on a set of 2D sections. The range of structural properties of such heterogeneous media reconstructions can be very broad, so that different classes of heterostructures can be simulated by this method, for example, Voronoi tesselations (mosaics). In the present state of random heterostructures physics one can formulate a large number of unanswered fundamental questions about the stochastic reconstruction approaches with specified correlation functions such as: What does the 2D analyses tell about the actual 3D structural details? What is the minimal amount of experimental data (microphotographs of 2D sections) necessary to reconstruct the actual medium, namely the porosity and the correlation functions?, How to quantify the nonuniqueness of the reconstructed systems? The solutions to these illustrative problems all require explicit quantitative knowledge of the geometry and topology of the fine structural (local) details of the material medium. Computing the effective permittivity of random heterostructures from first principles is still one of the most challenging problems in materials science and is currently one of the limiting factors in ab initio predictions of the dielectric behavior of composite materials. The use of case-by-case calculations with no clear heuristic guidelines seems misguided in the face of this simple phenomenology. Such numerical approaches may ultimately lead to a systematic means for a multiscale analysis of random dielectric heterostructures. The resulting physics is remarkably complex and multifaceted. This leaves us with a richness of phenomena which is in fact necessary to confront the unexpected observations about electromagnetic transport in mixed media.

3. Historical perspective The idea that one can create new materials by mixing given ones on a macroscopic scale is certainly not new. Before plunging into the details of our numerical approach, it may be worthwhile to pose a general question: Does one need to know anything about the history of a research area such the dielectric behavior of heterogeneous media in order to appreciate the subject matter? Most of us are complacent about quoting the usual sloppy misattributions of published work even if we are finicky about the details of the experimental results and theoretical proofs [64]. We note in passing that aside from the human interest involved in biographical studies, there may be some intellectual value in retracing the way physical ideas have developed. This development is often messy, however. Occasionally good ideas emerge prematurely in obscure places and are forgotten for a time, only to be rediscovered independently. But, in the end, one is often curious to know where the currently accepted ideas came from.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

383

In this section, we review some important milestones concerning our understanding of the dielectric properties of heterogeneous materials. The science of dielectric heterostuctures has a long and distinguished history with much engineering literature devoted how to deal with these materials. In the past, many models have been discussed in a number of approximations in attempts to understand the underlying physics governing this issue. The historical development of research work on the macroscopic behavior of heterogeneous media has been documented in several reviews. A summary of contributions made by the end of the 1970s was presented by Landauer [1], and a more recent review article has been presented by Bergman and Stroud [3] in 1992. If the reader is interested to learn more about homogeneization formalisms, Ref. [65] is a good place to start; there you will also find additional references for much of the work described in this article. It is also worth keeping in mind that the vast literature related to the subject cannot be mentioned in a single review article. This paper does not attempt to be a comprehensive review of the field of dielectric heterostructures, even as regards its computational aspects. We used our own judgment and prejudices to select the subjects that we considered more relevant, and apologize beforehand for any omission. Only the basic results will be presented; the interested reader may then find the relevant details in the published literature. Formulas for the permittivity of heterostructures were originally discussed by Maxwell and Faraday in the middle of the nineteenth century as part of their comprehensive studies of electricity and magnetism. Faraday (1857) in a classical study made the first attempt at exploring the optical effects originating from the dispersion of fine metal particles in some host media [66], and Maxwell (1873) investigated the simplest kind of two-component composites consisting of spherical inclusions of material of type one embedded in a medium of material of type two [67]. By considering a single sphere in a continuum, Maxwell solved the Laplace equation for the potential both inside and outside the sphere, and using the principle of continuity, he obtained an expression describing the variation of the potential in the continuous phase due to the presence of the sphere. Subsequently, he obtained the effective conductivity of a dilute dispersion of spheres by regarding the dispersion itself as a sphere, and setting the effect of the dispersion on the potential in the continuous phase outside of the spherical boundary to that obtained for the case in which the sphere consisted of a single inclusion only [67]. These approaches for studying the dielectric behavior of matter is to adopt continuum electrodynamics, the macroscopic view of electromagnetic fields [68–70]. Not longer after those pioneering works, early researchers found approximation formulae by making ad hoc guesses about the effective electromagnetic properties from the properties of their constituents, i.e. mixing laws, provided that the different constituents of the composite do not react with each other under the mixing condition and with only small diferences between their respective dielectric permittivity. Briefly stated, within mixing law approaches the effective permittivity of the inhomogeneous medium is conveniently equal to a volume average of the local permittivity. The relative merits of each formula have been the subject of much discussion, but none is outstandingly successful. Modelers are often forced to rely on mutually inconsistent results drawn

384

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

from disparate sources, and even employing intuition and analogy to guess at effective permittivity values. Many papers reported impressive agreement between these theoretical conjectures and their experimental results [1,2,10]. However, people mostly examined the very dilute extreme cases with respect to inclusions. Moreover mixing rules have been criticized from the philosophical point of view and tthere are apparently conflicting claims about the accuracy and the range of applicability of the various formulae. Our [13–20] and other [3,26,29] previous studies have demonstrated that mixing laws can be used to fit experimental data, but cannot be generalized to make accurate predictions. In fact, this is not surprising to consider that the effective permittivity is not uniquely determined by the volume fraction and the permittivity of the different constituents of the composite material [71,72]. However, proper determination of the form of the particle–particle and particle–matrix interaction potentials is necessary to compute the ’’multipole expansion’’, which isolates the physics much more robustly [46]. There are many other notable names such as Mossotti (1850) [73], Lorenz (1880) [74], Lord Rayleigh (1892) [75] who investigated the conductivity of periodic assemblies of conducting inclusions, i.e. parallel cylinders in a square array and spheres in a regular rectangular array, of uniform size and of a given conductivity packed in a conducting host of a different conductivity, Clausius (1894) [76], Maxwell Garnett (1904) [77], who discovered the eponymous formalism for isotropic–indielectric mixtures, Wiener (1912) [78], Lorentz (1916) [79], Fricke (1924) [80], Lichteneker (1926) [81], and Bruggeman (1935) [82] who published a seminal paper in which the notion of effective-medium approximation (EMA) was popularized for the first time. Both, Maxwell Garnett [77] and Bruggeman [82] analyses become exact for a composite material composed of a small volume fraction (dilute limit) of particles randomly dispersed into a host matrix, but the latter shows eventually new features, for example, metal–insulator transition if the particles are conducting and the host is a dielectric; in those materials the transition takes place via an inhomogeneous regime in which percolation effects determine the transport properties. These models provide an understanding of average properties and limiting behavior, but are relatively inflexible. Although the seeds of the idea of EMA also appeared in 1806, in a classic paper by Biot and Arago [83], the rigorous EMA formalism had to await the advent of the Maxwell equations. The virtue of the methods based on the EMA is that they are not limited to low volume fractions of inhomogeneities. In fact these self-averaging (decoupling) procedures underly the mean-field character of the EMA. While the EMA has been successful in explaining the overall trend of the volume fraction dependence of the effective permittivity observed in matrix–particle composites, it must be stressed that the results are somewhat ambiguous in view of the various approximations involved in coping with such complex systems. First of all, the analysis is one of continuum physics: EMA approaches assume that each constituent is surrounded by the same effective medium. It assumes that the local electric and magnetic fields are the same in the volume occupied by each component in the composite. The analysis is done in the approximation of noninteracting inclusions (each inclusion is subject to the same mean-field, unperturbed by the presence of other inclusions). The accuracy and range of validity of EMA is not easy

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

385

to establish. In principle, this approximation is rigorous at small volume fraction of inclusions (provided that their mutual positions are uncorrelated so that cluster-type configurations are excluded). Thus, EMA analyses are insensitive to the detailed structure of the composite material, for example, connectedness and clustering of one phase for two-component random heterostructures. In other words, in the effective medium the energy density is homogeneous by construction. Another problem of principle with this philosophy is that, although it produces the correct compositional dependence for the permittivity at low volume fractions of inclusions (dilute concentration limit), EMT as previously formulated certainly break down at higher volume fractions of inclusions since they are based on a dipolar analysis, and do not take into account the multipole correlations, i.e. many-body problem inherent in the electromagnetics of the particular system. It is also worth noting that what the matrix is and what an inclusion is are determined by the physical situation of the specific problem at hand, and it cannot be arbitrarily interchanged just to obtain an agreement with some measurements. It is also worth noting that Stroud and Pan [40] have extended the EMA to finite frequencies, using a full multipole expansion to treat scattering from particles. Despite their undoubted value, one must be extremely cautious in drawing conclusions concerning the dielectric behavior of matrixparticles composites from investigation of macroscopic effective permittivity by means of EMT arguments alone. Exact calculations on regular arrays of conducting spheres have shown that the effective permittivity becomes infinitely great as the spheres approach contact [21,47]. One lesson learned from EMT is that even if the geometry is known, the problem of predicting the dielectric behavior of simple twocomponent heterostructure models is difficult. While the correlation length
usually characterizes the range of validity for EMT, this is a poor approximation in the vicinity of a scattering resonance [3,84]. When such resonance processes occur, it is necesssary to formulate the calculation of the effective permittivity by separating out the resonance state and treating it explicitely [3,38]. On a similar perspective, the Bergman–Milton [36,85] spectral formalism permits to separate the effect of the dielectric properties of the constituent materials from the effect of the geometry of a two-component composite. Appropriately, Bergman and Milton appreciated the fact that the effective properties are given by an integral transform of a function, known as the spectral function (a special type of Stieljes function described by a series of simple poles and residues), that depends only on the geometry of the composite and are independent of the specific material properties of the inclusions, i.e. once it is known, it provides a solution for any type of inclusion that has the same geometry [86,87]. Later, the development of the generalized EMT by McLachlan [88,89] has been also a strong impetus for research in the dielectric properties of inhomogeneous materials. McLachlan’s model, developed heuristically, has its own specific peculiar features and offers practical advantages in predicting composition dependence of the effective permittivity of composite media. Although McLachlan’s formulation has not been derived by mathematical arguments, it is reasonable because it reduces to Bruggeman’s EMT and has the mathematical form of the percolation model in certain limits. Even after 20 years of active research, there is still controversy as to the

386

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

theoretical justification and experimental verification of this model. While several experimental studies [2,88] showed in some cases that McLachlan’s model does reproduce certain gross features of the qualitative composition dependence of the effective permittivity of real compounds, it is generally believed that the most severe approximation of this model is the conspicuous absence of first-principles and selfconsistent arguments of its derivation. More specifically, the experimental results raise a wide variety of questions such as whether this model is relevant to any type of structural disorder, whether the rapid rise in permittivity up to a singularity is due to a percolation effect, and what is the universality class, if any. Very recently, Brosseau [90] has argued, from measurements of the complex effective permittivity, that the microwave relaxation in carbon black filled polymers can be accounted for by incorporating Jonscher’s model in the McLachlan method. That the effective permittivity of these materials is sensitive to the percolative process of carbon black particles is a reasonable idea and in fact expected. To incorporate the information on percolation in the phenomenological approach of McLachlan, it requires the knowledge of the critical exponents s and t. Other attempts toward improving EMT introduce a set of phenomenological parameters, for example, the coated spheres technique [84,91,92]. However these phenomenological models are required to invoke additional (often arbitrary) assumptions to reduce the huge number of physical and microstructural parameters, i.e. the distributional statistics of the constituent phases, to a manageable few. Quite often, we do not have the detailed knowledge of this information, but only the values of the permittivity of the constituents of the composite, the corresponding volume fraction, and in some special cases symmetry properties, for example, rotational invariance of the permittivity tensor for isotropic composite materials. Faced with this situation, many researchers turn to completely different approaches in their quest to evaluate the effective permittivity of composite materials and focus their effort in calculating bounds on the allowed values of the permittivity. A simple but very appealing picture for analytical bounds of the effective permittivity of a binary lossless composite in terms of its constituents parts goes back to Wiener [78]. Subsequent theoretical progress was made by several investigators who attempted to optimize these bounds by introducing ‘‘details’’ of disorder as well as certain symmetries, notably Beran [93], Hashin and Shtrikman [94], and Bergman [36]. While Beran [93] and Hashin and Shtrikman [94] make use of a classical variational principle to dertermine bounds for the effective permittivity of macroscopically isotropic random media, Bergman exploits the analytic properties of the effective properties as a function of the constituent properties. Some practical applications of these bounds have been proposed [61]. For completeness, it would be appropriate to mention that rigorous bounds on the complex effective permittivity for lossy two-component composite media have been derived by several investigators including Schulgasser and Hashin [95], Milton [96,97], Sawicz and Golden [98], and Bergman and Stroud [3]. These bounds are represented by limited domains in the complex plane inside which the allowed values of the permittivity exist. ‘‘Improved bounds’’ depending non trivially upon correlation functions, i.e. containing information beyond that embodied in the volume fraction of the different constituents, have been proposed by several investigators.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

387

From the experimental standpoint, it is generally a very laborious task to determine the correlation functions which characterize the structure of the composite material, for example, Corson [99] has made measurements of three-point correlation functions for Pb–Al and Pb–Fe mixtures. An earlier discussion [6] has demonstrated that the bounds become progressively narrower as more microstructural information is included, and may serve as a guide in identifying appropriate statistical descriptors. A new twist to the story appears in 1957 with the advent of the percolation concept introduced by Broadbent and Hammersley [100]. Later, a number of workers transformed Broadbent and Hammersley’s qualitative explanations in the context of critical phenomena which has been known as the percolation theory [101–103]. Historically, the percolation model was proposed as a conceptual aid to describe transport properties of disordered systems with a large disorder. Both the formalism and many of the attendant ideas related to percolation lie at the core of a significant fraction of contemporary research in the study of random systems. In physics, percolation has provided a great deal of insight in the study of microscopic [104] as well as macroscopic phenomena [105]. This concept finds application in a number of scientific disciplines and in engineering, for example, polymer gelation, transport in porous media, epidemic spreading, and occupies a prominent place in contemporary thinking about complex disordered condensed matter systems. The physics of the percolation transition itself has stimulated broad interest [101–103]. Let us first briefly review the percolation transport properties of the discrete lattice percolation model, in which the bonds (or sites) of the lattice are randomly occupied and all occupied bonds are identical. The key feature of the model is as follows: statistically speaking, it might be argued that above a critical bond concentration, there is an infinite cluster of connected sites, and below it there is not. Here it must be emphasized that we define the percolation threshold in a geometric sense, referred to in the literature as ‘‘connectedness percolation’’. As is by now commonly agreed the percolation model can provide a means to understand the nature and threshold of the insulator–conductor transition at a critical volume fraction of conducting inclusions into a dielectric matrix, for example, conduction in bulk carbon black-filled polymer with relatively high volume fractions of carbon black has been shown to be percolative [106–115]. But it should be added that the relationship between formation of a chemically or physically connected infinite cluster and appearance of conductivity threshold, is complex and dependent on details of the microscopic filler particles/ polymer chains interactions. A number of numerical simulations on 2D and 3D networks, and also in continuum, have revealed the underlying features of the mechanism for the connectivitiy transition in these systems and have been particularly insightful in elucidating aspects of the conduction threshold in composite materials [101]. Another characteristic property of these percolating systems is that close to a critical point, i.e. the percolation threshold, the dielectric behavior is not determined by the type of material but by its critical properties, which have universal descriptors, i.e. the universal critical exponents s and t. Near percolation threshold fc tiny modifications in volume fraction f of inclusions can induce large permittivity changes, as represented by power-laws "0  j f  fc js for f > fc and f< fc,

388

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

and "00  j f  fc jt for f> fc, for the real and imaginary parts of the complex permittivity, respectively. To the surprise of the researchers, the critical exponents were ‘‘universal’’, so that widely different systems had precisely the same exponents s and t. These values turned out to depend only on the coarsest physical properties, such as the spatial dimension, for example, s=0.73 and t=2 in 3D [103]. Another element of the scaling description is that the correlation length
diverges near the percolation threshold following a power law
 j f  fc j as f ! fc. The diversity of critical systems that can be described by universal relations indicates that experimental measurements on one system should yield complementary information as on another. Given the wealth of these results, some of the recurrent questions concerning the conduction threshold in composite materials can be addressed thanks to the percolation model. In particular, it became esssential to establish those general features of experimental and numerical results that are model independent to understand the physical nature of the conduction threshold in composite materials. Within this context scaling and universality concepts that arise in many circumstances in our physical world [103] appear particularly attractive for the understanding of the dielectric behavior of matrix–particles composites. The core physical idea of the scaling invariance is that, as the physical system moves toward a phase transition, it becomes increasingly dominated by large-scale fluctuations. At the critical point, the correlation length, i.e. the length scale of the fluctuations, becomes infinite and the system exhibits scale invariance. In spite over thirty years of theoretical research, however, we still do not have a complete understanding of the model itself, particularly finite-size and long-range correlation effects. Additionally, one should note that EMT of disordered solids are likely to be poor approximations near conduction threshold [3,13–19,90,109], for example, in the frequency range corresponding to the plasmon resonances in metal grains for metal-dielectric composite materials. Between the time of early experiments and subsequent analytic approaches for modeling the effective permittivity of heterostructures this field lay dormant. The reawakening of the field came as a consequence of the availability of powerful computers and the improvements in the numerical models, and it has allowed us in recent years to understand and predict the properties of some materials systems from first principles. The situation is particularly true for multicomponent heterostructures. Computational results can contribute to the acceptance and rejection of EMT, and can also indicate directions in which new approaches should be developed. Much progress has been in this direction, with numerical calculations routinely being viewed as impartial referees that may eventually select the proper analytic description of a given model. Actually, the cross-fertilization between computational and analytical work in the area of effective properties of heterostructures is quickly growing. Over the past several decades, a number of theoretical and algorithmic techniques, such as the use of Fourier expansion [116], FDTD [27,28,31,41–43], fast multipole method (FMM) [44], Harrington’s method of moments [117], finite element [118–120], finite difference (finite integration theory) [121,122], variational formulation [123], and BIE [13–20,23,124] are among the techniques which have been vigorously pursued in computational electromagnetics. Details of the methods

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

389

may be found in the references given. At their most basic the full electromagnetic treatment of these methods involves the resolution of Maxwell equations which are coupled with specific boundary conditions. These new approaches make a number of specific predictions for the effective properties of matrix-particle composites. This kind of approach has not been free from ambiguities either. These arises from shortcuts imposed by the computational cost: modest size of the model systems, interface effects between the particles and the host matrix, etc. It is not surprising that no definitive numerical approach has been firmy established, given the complexity of the problem and the difficulty of isolating and identifying the physical parameters at work. Current ab initio methods allow for determination of the effective permittivity of matrix-particles composites. As was mentioned earlier we have developed, in a recent series of papers [13–20], a mathematical model for determining the solution of boundary-value problems defined on domains with specified geometries in periodic composites. The formal development of this model leans heavily on the use of Green’s theorem and boundary integral equations are solved for the local field with appropriate periodicity conditions on a unit cell of the structures using finite element. Finite element method modeling is a numerical technique in which the system to be analyzed is subdivided into a number of elements with finite dimensions. These elements are interconnected by means of nodes representing points in the macroscopic system, which can be used to predict macroscopic behavior. Technical details about calculations related to such BIE approach will be presented in the following section. Applying Green’s theorem, we can obtain integral representations of electromagnetic fields from Maxwell’s equations. The idea of deriving boundary-integral equations from spatial Green’s functions is absolutely classical in electrostatics [69,70]. Our BIE approach to dielectric composite materials is similar in some respect to that employed by Ghosh and Azimi [23] in a related context, for example, case of nonabsorbing particles in a lossless matrix. We note that this treatment has been also employed by Sihvola and coworkers for related problems [27–29]. Numerous advantages of the BIE approach have been shown through different applications [13–20]. Later in the paper we shall show that integral methods offer several advantages, as compared to other numerical techniques. For instance, only the region where the sources are located must be discretized. This allows to avoid the meshing of the air, and the related need of dealing with regularity conditions for unbounded domains. Despite these important aspects, integral methods are not very popular, since they require the generation, the storage and the inversion of dense matrices [124]. In addition, this analysis do not cover mixtures of inclusions of different shapes (polydispersity). During the past two decades, there has been a huge experimental and theoretical efforts devoted to the study of the electromagnetic properties of random materials: Understanding the disorder effects arising from quenched randomness is one of the central problems in condensed matter systems. The two main families of computational models employed are off-lattice models and lattice models. The former are typically more versatile and realistic while the latter are usually more computationally efficient. The heavy computational demands of ab initio permittivity calculations restrict their use to the study of small condensed matter systems and relatively few configurations. More-

390

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

over, the treatment of this problem using first principles is increasingly cumbersome, if not impossible in many situations because of the complexity of the interparticle interactions. The understanding of the physics in turn allows for such calculations to be performed using the BIE approach albeit in an approximate manner, with the obvious requirement that the dimensions of the cell which must be larger than the correlation length
. Despite this recent focus, the computational evaluation of the effective permittivity in the presence of strong correlations (say, due to the long-range multipole interactions) is still rather poorly understood, even though some important attempts have been made [6,125]. There has been also a significant advance in recent years in the investigation of finite-size effects. The downside of these approaches is that, in order to obtain converged simulation results for a permittivity value of random heterogeneous materials, one needs to employ large system sizes, which makes this approach computationally very demanding. Since numerous systems of industrial and technological interest fall precisely in this category, it is important to develop alternative methods for calculating the n-point correlation functions, which would improve the accuracy of the BIE formalism, while retaining its computational simplicity. These types of computations facilitate the development of physical insight into the relative importance of the physical parameters in the system. In this short historical survey, we have spanned the exciting two and a half-century period from Faraday’s memoir of 1857 to the most recent research reports dealing with computational electromagnetics. Despite the extensive work, a number of big questions about the dielectric behavior of heterostructures which were raised by the original discoveries of Faraday, Maxwell Garnett and Rayleigh, to cite but a few, remain. Even though many of the experimental data have been available for a long time, one cannot help but notice some lack of perspective in the early interpretations, characterized by the emphasis on only certain of the physical mechanisms. This is not to say that we have not learned an enormous amount about both from experimental and theoretical studies : but to isolate what is crucial to the physics from its irrelevant clothing is the central task. In order to improve the understanding of the physics of dielectric heterostructures a much broader point of view is necessary, one that is able to put together the whole body of experimental data, in addition to the modern concepts of condensed matter physics.

4. Numerical procedure The purpose of this section is to recapitulate the salient features of the formal rigorous BIE method employed in the present investigation for calculating the effective permittivity of matrix–particle composites. Along the way we pick up notation and acquaint the reader with the kind of thinking that underpins the principle of the numerical approach we will use. Some of the key equations are reviewed and this is followed by some computational details relevant for the test cases to be studied in the next section. Observe at the outset that a number of different effective transport coefficients in heterogeneous materials can be calculated using the same

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

391

computational method, for example, elastic moduli, fluid permeability, thermal conductivity, for reasons of mathematical analogy [126]. 4.1. Problem statement To establish notation: For the sake of convenience we shall assume that the composite medium, which is static in time, consists of lossless homogeneous inclusions of permittivity "1 embedded in a lossless host matrix of permittivity "2 . Since we are considering lossless dielectrics, the permittivity of the effective medium has to be real due to energy conservation [70]. For simplicity, the inclusions are assumed to be connected regions with smooth boundaries. A wide variety of inclusions can be considered regularly arranged on a 2D, or a 3D array. Figs. 1 and 2 provide the reader a schematic overview of the two regions associated with 2D (Fig. 1), or 3D (Fig. 2) lattices. We define a unit cell of the structure as a parallelepipedic (or cubic) volume which can generate the whole composite structure by applying appropriate translations. It is worth observing that the way of defining the unit cell is not unique, but previous work has shown that the different choices give equivalent results for our purpose [13–19]. The method is general for composites containing any number of component materials as long as the composite has a periodic lattice. In fact, symmetry permits us to describe the periodic composites by 18th of each unit cell. The BIE treatment is based upon the potential distribution in this unit cell. So far, no approximations have been made. In considering electromagnetic propagation in such mixed media, it is important to bear in mind that the validity of the QA implies an upper bound for the frequency of the electromagnetic wave probing the system. According to the condition l
and restricting the discussion for nonmagnetic (=1) and lossless (" real and positive) materials, we can evaluate, for heuristic purpose, the range of applicability of the QA with the following inequality:

Fig. 1. The finite element model used in the computation of the effective complex permittivity of twodimensional periodic composite materials. Boundary conditions are set on the limits of the unit cell structure. The applied field is oriented along the y-axis. The dielectric constant of the remaining space is "2 .

392

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

!  !m ¼

1 2c 2c 12 ð"Þ2 ¼ ð"Þ



c being the velocity of light in vacuum. To be more explicit, let us consider two examples illustrating this point: using a value of
=10 nm and "=3, one finds m ¼ !2m =1017 Hz, the inequality holds for frequencies up to the visible range of the spectrum and if
=1 mm and "=3, we get m=1012 Hz, as a result QA is a justified approximation up to the microwave frequency range. In addition, we assume the inclusion size to be small compared to l so that there is no skin effect. The physical implication of this results is that we can achieve a macroscopic description of the composite as a homogeneous local continuum. For a periodic array, wavelength should be much longer than both the inclusion size and the lattice spacing. Thus, we treat the material medium as a homogeneous equivalent material. The effective (relative) permittivity " can be defined as E 2" ¼

1 O

ð

"^ jrVj2 dO O

Fig. 2. Boundary conditions related to the 3D configurations investigated in the numerical computation. The dielectric constant of the remaining space is "2 : (a) isolated particle of permittivity "1 , (b) fused particle of permittivity "1 .

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

393

ð 1 EdO denotes the average field which O O depends on the applied potentials, V is the potential, E is the local electric field in the cell and O is the volume of the medium. Note that this definition ensures that if the composite medium were replaced by a homogeneous medium, subject to the same boundary conditions, the rate of dissipation would be the same. Consequently, higher electric multipole interactions are treated exactly. This is of fundamental importance because higher multipole interactions become important when the inclusions approach contact. This method goes beyond the mean field approximation which is at the basis of many effective medium laws. Having made the assumptions as indicated above, we are able to solve the problem at hand exactly within BIE theory. Briefly stated, we calculate the local potential in response to an external applied electric field and determine the effective permittivity of the composite from the calculated potential. The effective permittivity is defined as the permittivity of an isotropic homogeneous material which reproduces the complex capacitance when inserted in place of a sample of the composite material as the filler of a parallel-plate condenser, across which an oscillating electric field has been applied [3,13–19,23]. It is quite important to remark that we also assume that the dimensions of the condenser are much larger than
but sufficiently small compared with the wavelength l. We use Green’s functions technique to derive integral equations for these structures and to transform the problem to a form in which efficient boundary element discretization can be used and a matrix inversion problem to be solved. To proceed with the approach, consider the boundary-value problem which is solved by using a six-step procedure as follows: (1) construction of the surface meshing of the unit cell, (2) generation of the inclusion with specification of volume fraction, size, aspect ratio and orientation, (3) application of prescribed boundary conditions given in terms of the potential (Dirichlet) or the electric field (Neumann), (4) specification of permittivity of the constituents, (5) calculation of the electric flux through each surface of the unit cell, and (6) calculation of the effective permittivity for a given direction of the applied electric field. It is convenient to display the relative permittivity " compared to the vacuum permittivity "0=8.85 1012 Fm1. The absolute permittivity, if desired, is obtaining by multiplying the relative permittivity by "0.

where "^ is the local permittivity value, E ¼

4.2. Algorithm We now specify the procedure to calculate the effective permittivity. Once a continuum quasistatic model is adopted, the main task consists in designing and implementing the model into a reasonably efficient computational scheme for finiteelement simulations. A valuable starting point for attacking such a problem is through a first principles analysis in electrostatics, namely by using Laplace’s equation, i.e. !.("!V)=0, where V is a potential distribution inside a spatial domain  bounded by a closed surface S, and which has no free charges or currents. Upon using Green’s theorem [46,69], we can write the local potential V(M2O) in terms of

394

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

V(P) and of the normal derivative @V@nðPÞ with P being any point on the boundary S (with no overhangs) of :  ð  4 @G @VðPÞ VðMÞ ¼  VðPÞ G ds ð1Þ A S @n @n where A stands for the solid angle under which the point M sees the oriented surface S, n is the normal unit vector oriented outward to S, ds is a surface element of S and G denotes the Green’s function. Refering again to Fig. 2a, we have represented an arbitrarily shaped homogeneous inclusion occupying a volume O1 and permittivity "1, which is embedded in a homogeneous matrix of volume O2 and permittivity "2 . At this point it is worth recalling that absence of electric charge density and of conduction current density are tacitly assumed through the analysis. Given these assumptions, Eq. (1) leads to    ð  4 @G @V G V¼ V ds; ð2Þ A S1 @n @n 1 for domain 1, and    ð  4 @G @V G V¼ V ds; A S2 @n @n 2

ð3Þ

for domain 2. We must also specify boundary conditions for the potential across boundaries: the potential is subject to the boundary condition     @V @V "1 ¼ "2 ; ð4Þ @n 1 @n 2 by virtue of the conservation of the normal component of the electric displacement vector at the interface, provided no double-layer exists at the interface. Consequently, we have to solve the above two integral Eqs. (2) and (3) to evaluate numerically the electrostatic potential distribution. For that purpose, the implementation of the boundary-integral element method [124] consists in dividing the boundaries into finite element and for each finite element, the calculation is carried out by interpolation of V and @V @n with the corresponding nodal values X lj Vj V¼ j

and   @V X @V ¼ lj @n @n j j where lj denote interpolating functions (see Appendix A for details) [123,124]. Consequently only the surface bounding the volume is meshed. Following this way, integral equations are transformed in a matrix equation which is numerically solved

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

395

using the boundary conditions on each side of the unit cell as displayed in Fig. 2a. To solve these matrix equations we use a matrix inverse iteration approach based on the Galerkin method of weighted residuals. One salient feature of our method is that only the boundaries of the geometry need be discretized, not the entire volume, which has for effect to reduce the memory space required for manipulation of data, but the matrix equation to solve is asymmetric and full. It is worth emphasizing that the entire procedure outlined above, including the creation of the mesh, makes use of the commercial software package PHI3D [127]. It is worth observing that the fully automatic generation of highly quality meshes plays a crucial role in the finite element analysis, influencing both ease of use and accuracy of a software package. The field calculation package PHI3D permits the closely controlled generation of finite-element meshes through the use of input files containing complete instructions for node-by-node and element-by-element mesh specifications, along with imposition of boundary conditions. One basic question is as follows: How the number of nodes in the mesh can affect the quality of the obtained results? We tested coarse meshes with meshes of increasing fineness and found that above a mesh size of a few thousands elements the numerical results for the permittivity were constant and reproductible. For the numerical solution we discretized the configuration investigated here with a number of elements ranging between 5000 and 10 000, that means a linear system of equations with typically 20 000 unknowns has to be solved. The total polarization characteristics of the configuration of inclusions are then available by matrix algebra. Some comments are in order regarding this calculation scheme. In Fig. 2a, the inclusion is isolated (low concentration) and thus the medium of permittivity "1 cannot intercept the sides of the parallelipipedic cell. In that case the effective permittivity, in the direction corresponding to the applied field, is calculated using the following relation ð

  @V V2  V1 S; "2 ¼" @n e S 2

ð5aÞ

where V2V1 denotes the difference of potential imposed in the z-direction, e is the composite thickness in the same direction and S is the surface of the unit cell perpendicular to the applied field. If the inclusions are allowed to overlap with other (high concentration regime), the region of permittivity "1 can intercept the sides of the parallelipipedic cell (Fig. 2b). In that case, we have to take into account the electric displacement flux through the area S1 for the calculation of the effective permittivity in the direction corresponding to the applied field. Then, Eq. (5a) needs to be changed in ð "2 S2

  ð   @V @V V2  V1 þ "1 ¼" ðS1 þ S2 Þ: @n 2 S1 @n 1 e

ð5bÞ

396

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

Finally, it should be noted that the permittivity was calculated for a 1V potential difference imposed in the z-direction: constant potentials V1=0 V on the z=0 face of the parallelipiped and V2=1V on the opposite face (Fig. 2a,b). We note parenthetically that, in the above development, we have tacitly used the fact that the permittivity of the two components of the composite medium is a real number. It is worth emphasizing that the generalization of this calculation scheme to complex permittivity (represented by the complex number "="0 j"00 , where the real part, "0 , is associated with the ability of the material to store electric field energy and the imaginary part, "00 , reflects the ability of the material to dissipate electric energy in the form of heat) of composite systems, i.e. lossy materials, is straightforward and has been given in Refs [14,16–19]. Generally the permittivity is a tensor quantity. For anisotropic systems, for example, considered in the next section, the tensor components are resolved by applying the external electric field along the principal axes of the sample. As mentioned earlier, it should be noted that the BIE method gives an accurate description of the electric potential by taking into account edge and proximity effects even for a high concentration of inhomogeneities. Therefore, this numerical value does not suffer from the disadvantages of the traditional boundary-value approach. Although the fully 3D problems are in a sense the most fundamental ones, the reader will have no difficulty formulating analogous composite structures in the context of 2D models [14]. Another important and generic characteristic of the BIE method is that only the boundaries of the geometry need be discretized which has for effect to reduce the memory space required for manipulation of data, but the matrix equation to solve is asymmetric and full. We also note that, in contrast to other computational schemes reported in the literature, the field calculation package PHI3D uses a semi-automatic meshing system. Finally we emphasize that the computational time to run our simulations on a HP model 712/80 workstation are linked with the details of the meshing of the geometry considered and ranges from a few minutes to a few hours for calculating the effective permittivity of a typical 3D configuration. The full appreciation of the magnitude of the computational savings afforded by the BIE method has been already discussed [13–15]. In the following section, we present 2D and 3D computation results that illustrate the utility of the ideas in the previous section. To assess the success of the simulation protocol we checked for consistency between the effective permittivity of composites obtained by this numerical approach and that predicted by several analytical methods. As will be shown later, the results obtained using the simulations agree well with the closed form expressions from diverse theories. This agreement inspires confidence that the simulations will correctly predicts results that are difficult to calculate in closed form using analytical theories.

5. Numerical examples and discussion In this section we systematically go through the dielectric properties of several relatively simple two-component periodic (2D and 3D Bravais lattices) heterostructures

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

397

using the numerical procedure just developed and explain how they are related to the nature, size, shape and orientation of the constituents. One goal of this extensive ab initio computational study will be to compare the method with previous work in several cases. Some of the results presented in this section have been briefly reported before. Here we present a complete exposition of these results, as well as many new results related to this problem. Finally we discuss a number of issues that must be addressed in order for a dielectric heterostructure to be manufactured as a practical device. 5.1. Two-dimensional simulations For conceptual purpose and to put our numerical results into perspective, it will be useful to begin by considering results of 2D simulations. This situation can model media whose material boundaries are cylindrical surfaces with generators parallel to one axis, for example, oriented fibers in a polymer matrix [110]. 5.1.1. Square lattice and 2-pole approximation As a check on our numerical code we have validated it for a simple 2D geometry for which analytical expressions can be derived [22]. Let us consider the situation where a 2D square lattice is formed by an array of cylinders, which may overlap. Fig. 3 shows the geometry of the unit cell of the square lattice with Cartesian spatial coordinates x, y and z. The structure consists of an infinite cylinder ("1 ; 1 ¼ 1) of radius r with generator parallel to the z axis embedded in a host matrix ("2 ; 2 ¼ 1), i.e. the intersections of the cylinders with xy plane form a periodic two-dimensional structure (simple square lattice, with e=1). This lattice is characterized by circular or piecewise circular boundaries between the materials. Note that this geometry

Fig. 3. Notation related to the configuration of the two-component 2D composite system (simple square lattice) investigated. The geometry is symmetrical in both the x- and y-directions. The unit cell is assumed to be a square of length e=1 containing an infinite cylinder of radius r along the z-direction and of permittivity "1 . The dielectric constant of the remaining space is "2 .

398

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

renders the two-component composite materials to be translationally and rotationally invariant. We denote by f the fractional occupancy of the constituent 1: f ¼ r2 below the percolation (0 4 r 4 12, non-overlapping cylinders) and fp¼ ffiffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  1  pthreshold 2 2 2 r  4r arccos 2r þ 4r  1 beyond the percolation threshold ( 12 4 r 4 22, overlapping cylinders). Note that in the latter case, we exclude the four sections of the component of permittivity "1 exceeding from the unit cell. To facilitate the comparison with the results of Ref. [22], we used the following set of parameters: "1 ¼ 1  3i and "2 ¼ 5  8i. In Fig. 4a,b, we show the comparison of the effective permittivity computed for two-pole approximation (Bergman–Milton

Fig. 4. (a) Real part of effective complex dielectric constant evaluated by the BIE method and the 2-pole Bergman–Milton approximation [36,85] as a function of the radius of the cylinder r. The dielectric constants of the two-component composite material are "1 ¼ 1  3i and "2 ¼ 5  8i, (b) Imaginary part of effective complex dielectric constant evaluated by the BIE method and the 2-pole Bergman–Milton approximation [36,85] as a function of the radius of the cylinder r. The dielectric constants of the twocomponent composite material are "1 ¼ 1  3i and "2 ¼ 5  8i.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

399

theory [36,85]) and the BIE method versus the radius of the cylinder. Figs. 4 and 5 are for the real and imaginary part of ", respectively. The solid dots are our numerical results and the crosses are the numerical results from the data reported by Liu and Shen [22]. The agreement of both sets of results is good. As can be recognized from these figures, the agreement of data on the entire range of r, indicates that our model is consistent with the Bergman–Milton theory [36,85]. 5.1.2. ‘‘Random’’ configurations Up to this point we have not paid any attention to composites which have some element of randomness. This brings us to the question concerning the relation of the spatial randomness to the dielectric behavior. The 2D lattice can serve as a guideline structure to demonstrate the modifications of permittivity that can be caused by ‘‘random’’ distributions of particles using the method of superposition of periodic structures [15]. The basic advantage of this approach lies in its ability to describe disordered heterostructures with terms (and accuracy) equivalent to those with which state-of-the-art methods address simple periodic structures. Like any

Fig. 5. (a) Geometry of the unit cell of the two-dimensional composite material investigated. The unit cell consists of a square diamond of permittivity "1 and half-diagonal inside a square of permittivity "2 and side 2. (b) Random distribution, (c) regular distribution, and (d) agglomerate distribution.

400

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

sampling method, the superposition of periodic structures approach is effective as long as the variance of the pertinent physical property is sufficiently smaller than the mean; otherwise a large number of structures needs to be included. In turn, this variance reflects the dependence of the property of the structure, for example, clustering, through its correlation functions. Moreover, for three-dimensional structures, the dimension of the matrix systems to be solved becomes very large. Thence the resulting CPU times to obtain solutions increase dramatically. In the case of 3D random composites the effective permittivity is computed by considering the equivalent periodic material (with identical inclusions oriented in the same direction) and taking a statistical mean of the permittivity in the three directions x, y and z:  1 " ¼ "x þ "y þ "z ð6Þ 3 This procedure may be justified by the fact that these media are macroscopically homogeneous and isotropic as previously mentioned. Here, the effective permittivity is a scalar parameter that can be derived from the effective permittivity tensor obtained for an anisotropic medium corresponding to periodic composites with oriented inclusions. In terminating this subsection it should be mentioned that while the BIE method is exact for periodic structures, its application to random media through Eq. (6) consitutes an approximation. It gets worse at high volume fractions of inclusions. The immediate focus here was to determine the influence of the topology of the clusters on the effective permittivity. It should be noted that in our experiments, the background relative permittivity is that of free space ("2 ¼ 1). Impenetrable inclusions were considered for simplicity. The materials being non-lossy, their permittivities are real numbers. To investigate the influence of the component random arrangement in space, several simulations were carried out. The unit cell of the 2D composite material is presented in Fig. 5a. This unit cell consists of a square diamond of permittivity "1 inside a square of permittivity "2 and side L=2. The 2 corresponding volume fraction of the inclusion phase is given by f ¼ 4 . To take into account the effect of the inhomogeneities distribution on the effective permittivity, the domain of permittivity "2 is 100 times duplicated and the inclusions of permittivity "1 are arranged randomly with concentration equals to C (see Fig. 5b). We define this concentration by the ratio of the number of inclusions to the total number of unit cells of permittivity "2 . Thence the resulting volume fraction of the random composite is expressed by f~ ¼ Cf. Examples of different kinds of distributions are displayed in Figs. 5b–d. The first pattern (b) describes a random arrangement of inhomogeneities (strong disorder). The second pattern (c) concerns inclusions that are regularly distributed in the host medium (low disorder). In this case, the arrangement is quasi-periodic. The last pattern (d) deals with agglomerate distribution of inhomogeneities. Inclusions aggregate to form isolate clusters. The effective permittivity is computed by the BIE method for these types of distributions and for two sets of permittivity ("1 ¼ 2, "2 ¼ 1) and ("1 ¼ 5, "2 ¼ 1). The resulting effective permittivities are compared with those obtained from periodically-arranged inclu-

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

401

sions, i.e. Fig. 5c. The effect of the distribution on the dielectric properties of the composites is observed in Figs. 6a,b. For a low contrast ratio between the permittivities of the background and inclusions, the effective permittivity of the composite is not affected by the inhomogeneities distribution. Increasing this contrast ratio shows the influence of randomness. The deviation between the results obtained for

Fig. 6. (a) Influence of the distribution of inhomogeneities on the effective permittivity. Inclusions with permittivity "1 ¼ 2 are embedded in the background matrix with permittivity "2 ¼ 1. The effective permittivity is computed by the FE method for the different types of distributions investigated. The open circles, solid squares, and the dashed line correspond to random, regular, and agglomerate distributions of inhomogeneities, respectively. The dashed curve describes the effective permittivity of the equivalent composite with a periodic arrangement of inclusions. (b) Same as in Fig. 6a for inclusions with permittivity "1 ¼ 5 which are embedded in the background matrix with permittivity "2 ¼ 1.

402

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

the three distributions investigated and those computed in the case of a periodic arrangement is large for a strong disorder and when inclusions form isolated clusters. While the results of such simulations are qualitative in nature, care must be used when the results are compared with other simulation data or experimental data. Can such numerical data reveal some information on the physical and physico-chemical mechanisms which are effective at the detailed quantitative level? Quantitative results may not be transferable between different types of randomness, but it may be useful to see which trends continue with more realistic systems. The results presented here indicate that the answer to this question is prematured. From this description one of the fundamental limitations to the study of random heterostructures becomes evident. Calculations with nontrivial (random) topology of inclusions embedded within a host matrix are extremely challenging, making it difficult to draw robust conclusions. We are trying to estimate the permittivity statistically from a finite number of samples, hence our estimates will be uncertain by an amount proportional to the square root of the number of samples [121]. More detailed numerical studies of realistic composites will be deferred to a subsequent publication. 5.2. Three-dimensional simulations Having assessed the accuracy of the theoretical method of calculating the effective permittivity of 2D systems, we now turn to the discussion of periodic composites with 3D lattices. The lattices that have been selected are the simple-cubic (sc) lattice, the body-centered-cubic (bcc) lattice, and the face-centered-cubic (fcc) lattice. 5.2.1. Nonabsorbing configurations 5.2.1.1. Spherical inclusions. Sphere represents an ideal starting point for study of more complex systems by virtue of its high symmetry. It is instructive first to consider a model with no absorptive processes. Hence, "1, "2 and " are real numbers. Consider equal-sized spheres fixed in a simple cubic array, a being the radius of the spheres. Fig. 7 shows a unit cell of the structure. For the purpose of simplicity, we assume, in the following, that all lengths ð‘; a; etc:Þ are dimensionless and that the side of the cell has the specific value ‘ ¼ 2 [13]. It is worth noting that if a  1 the particles act like isolated ones: they will experience only the external field and not the fields induced by the other particles. We call this case the isolated particle regime. The polarizability of the spherical particle can be deduced from the solution of the internal field of a spherical inclusion in a quasistatic field [46] 4 "1  "2  ¼ a3 "2 3 "1 þ 2"2

ð7Þ

Because of the periodicity of the structure and the spherical symmetry of the inclusions, the medium is isotropic. Therefore the depolarization dyadic L can be reduced to only one depolarization factor L characterizing all diagonal terms, i.e. L ¼ 13. By substituting this value into Eqs. (A1) and (A2) of Appendix A, it leads to Maxwell Garnett equation:

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

" ¼ "2

"1 þ 2"2 þ 2fð"1  "2 Þ "1 þ 2"2  fð"1  "2 Þ

403

ð8Þ

3

where f ¼ a6 is the volume fraction of inclusions in the mixture. It is a relatively simple exercise to show that Eq. (8) can be also written in the Clausius–Mossoti form (Appendix B): "  "2 "1  "2 ¼f ð9Þ " þ 2"2 "1 þ 2"2 What occurs when upon increasing the concentration f ? The distance of separation between two spheres decreases and particles will experience the local fields induced by other particles. The volume fraction f
¼ =6 ffi 0:523 corresponding to the limit of touching spheres (a=1) is the maximum packing threshold. Beyond that concentration, the geometry can be described according to Liu and Shen [22] excluding six segments of sphere from the unit cell. In this fused particle regime, the volume fraction is analytically calculated for a value of the radius a in the range pffiffiffi 1 4 a 4 2: f¼

a3   ða  1Þ2 ð2a  1Þ 4 6

ð10Þ

By taking into account the symmetry of the unit cell in the two regimes described earlier, the geometry is further reduced to one eighth of the microstructure for calculations by the BIE method. Our results are compared with those derived from Eq. (8) and with numerical data obtained by McPhedran et al. [21] and Tao et al.

Fig. 7. Notation relating to the unit cell of two-component periodic composite material investigated in the BIE model computation. The spherical inclusion, of radius a, with dielectric constant "1 is periodically arranged in a three-dimensional cubic lattice structure. The dielectric constant of the remaining space is "2 .

404

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

[116] for "1 ¼ 3 and "2 ¼ 1 (Fig. 8). In that case, it appears that the values of " computed by the BIE method agree satisfactorily with these previous calculations. One might argue that the agreement of the numerical data with Eq. (8), for f > f
, is somewhat fortuitous. This is confirmed by Fig. 9 which shows the limits of the Maxwell Garnett theory for volume fractions f higher than 0.4, when the ratio of the constituent permittivities is higher (we take "1 ¼ 30 and "2 ¼ 1). We also observe that our results are very close to those obtained by the Fourier expansion technique [116] over the entire range of volume fraction. The dependence of the composition on ", for f> 0.4, can be understood in terms of the multipolar character of the particle–particle interactions. As an interesting aside of these calculations, the results for the potential (resp. electric field) distribution of the lattice configuration considered for a volume fraction of spheroidal inclusions equal to 0.35 is shown in Fig. 10a (resp. Fig. 10b). 5.2.1.2. Ellipsoidal inclusions. Another interesting geometry to consider is to examine now an aligned ellipsoid composite that contains a periodic array of ellipsoidal inclusions in a host medium as shown schematically in Fig. 11. Composite ellipsoids consisted of nonabsorbing matrix containing a volume fraction f of nonabsorbing ellipsoidal inclusions. Fig. 12a and b show the details of the finite-element meshes. These meshes which amount to a mathematical description of the inclusion geometry, served as a foundation for our computational electromagnetic analysis. The unit cell can be described using the one obtained for spherical inclusions thanks to a tridimensional homothecy of ratio (‘; w; h), where ‘, w and h denote, respectively, the lengths of the sides of the parallelipipedic cell displayed in Fig. 11. As can be

Fig. 8. Volume fraction dependence of the effective permittivity of the three-dimensional periodic composite displayed in Fig. 7. Inclusions (permittivity "1 ¼ 3) are spherical and of volume fraction f in a matrix of permittivity "2 ¼ 1. The full circles are obtained by the BIE method. The solid curve corresponds to the Sillars equation [2].

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

405

seen in Fig. 11, the medium is anisotropic and the effective permittivity component "i in the direction i=x,y,z can be written as ni "i ¼ "2 þ ð11Þ ni 1  Li "2 where Li and i are the depolarization factor and the polarizability in the direction characterized by the index i. A standard result of electrostatics gives ð abc þ1 du pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Lx ¼ ð12Þ 2 2 0 ðu þ a Þ ðu þ a2 Þðu þ b2 Þðu þ c2 Þ where a, b and c denote the semiaxes of the ellipsoid. Note that Ly and Lz can be evaluated by interchanging b and a, and c and a respectively [46] . It is worth remarking that the depolarization factors do not depend on Pthe absolute values of a, b, and c, but only on their ratios. Moreover, the Lis satisfy i Li ¼ 1, i=x, y, z. If the scatterers are sufficiently distant from each other, their polarizability can be deduced from the solution of the internal field of a dielectric ellipsoid in a quasistatic field "2 i ¼ 0 ð"1  "2 Þ i ¼ x; y; z ð13Þ "2 þ ð"1  "2 ÞLi where 0 ¼ 43 abc is the volume of the ellipsoid. By substituting Eq. (13) into Eq. (11), the effective permittivity of the medium can be expressed as   ð"1  "2 Þf "i ¼ "2 1 þ i ¼ x; y; z ð14Þ "2 þ ð"1  "2 Þð1  f ÞLi

Fig. 9. Same as in Fig. 8 for "1 ¼ 30 and "2 ¼ 1.

406 C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

Fig. 10. Potential distribution for low (a) and high (b) concentrations of spherical inclusions in a simple cubic lattice configuration (Fig. 7).

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

407

Fig. 11. Notation relating to the unit cell of two-component periodic composite material investigated in the boundary element model computation. The ellipsoidal inclusion with dielectric constant "1 is periodically arranged in a three-dimensional cubic lattice structure. The dielectric constant of the remaining space is "2 .

Fig. 12. (a) Mesh used in finite element calculations of potential distribution for an ellipsoidal inclusion, (b) 18th of the domain is needed in this calculation due to symmetry considerations.

408

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

where f ¼ n0 ¼ 43 abc ‘wh is the volume fraction of the constituent 1. To simplify further the analysis, we first choose b ¼ c ¼ a2 (prolate spheroid) and we take the dimensions of the elementary cell as: ‘ ¼ 4 and w ¼ h ¼ 2. Fig. 13a and b show a comparison of the numerical results obtained by the BIE method for the effective permittivity of the medium in the x- and z-directions, respectively ("1 ¼ 30 and "2 ¼ 1) and those obtained from Eq. (14). They clearly show the limits of the analytical model for volume fractions f higher than 0.2 particularly; the effect being important in the x-direction for which the interaction effect is more stronger. Fig. 14a and b illustrate how the dependence of " versus f is affected by changing the semi-axis c of the ellipsoidal inclusion. As a model calculation, we take the following parameter values: ‘ ¼ 8, w ¼ h ¼ 2, "1 ¼ 30 and "2 ¼ 1. We observe from these figures that decreasing the semi-axis c by a factor of two has little effect on the permittivities "x and "z over the entire range of concentration. 5.2.1.3. Cylindrical inclusions. In a similar manner we consider the periodic composite displayed in Fig. 15. By taking into account the symmetries of the unit cell, we need only to evaluate the permittivity in the x- and z-direction since "y ¼ "x . The dielectric constant in the perpendicular direction to the cylinder axes is given by the Rayleigh’s formula

"x ¼ "2

"1 þ "2 þ fð"1  "2 Þ "2 þ "2  fð"1  "2 Þ

ð15Þ

  using Eq. (14) for prolate spheroid with Lx ¼ Ly ¼ 12   and Lz ¼ 2, where   1, van Beek [72] derived a general expression for the effective permittivity in the direction which is parallel to the cylinder axes "z ¼ "2 þ

1 ð"1  "2 Þð5"a þ "1 Þ 3 "a þ "2

ð16Þ

where "a denotes the apparent permittivity of the medium, i.e. ‘‘seen’’ outside by an inclusion. Its value differs from the permittivity of the host medium but must lie in the range "2 4 "a 4 ". To make our simulation concrete, the geometry of the unit cell characterized by a single parameter a by taking: H=8a, D=a, ‘ ¼ w ¼ 1, and h=8 (see Fig. 15). The 3 corresponding volume fraction of the constituent 1 is given by f ¼ a4 . Fig. 16a and d show data computed with the BIE method. At this point, a number of comments are in order. It is first interesting to observe that numerical data concerning "x are well represented by Eq. (15) for inclusions of low permittivity ("1=3), while for high permittivity ("1=30), there is a significant departure from Eq. (15) at high concentration levels. As concerns "z we note that numerical data for inclusions of low permittivity ("1=1) are well described by Eq. (16), taking an apparent permittivity identical to the matrix permittivity, while for higher permittivity ("1=30), Eq. (16) is

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

409

Fig. 13. (a) Volume fraction dependence of the effective permittivity in the x-direction of the 3D periodic composite displayed in Fig. 11. Inclusions (permittivity "1=30) are ellipsoids (b=c=a/2) and of volume fraction f in a matrix of permittivity "2=1. The full circles are obtained by the BIE method. The solid curves corresponds to Eq. (14) with a depolarization factor Lx=0.1736 [Eq. (12)]. (b) Same as in Fig. 13a for the effective permittivity in the z-direction, with a depolarization factor Lz=0.4132.

410

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

unable to correctly evaluate "z as shown in Fig. 16d over the range of concentrations explored. Besides the examples described in this work, other types of geometry of inclusions, for example, cube, disc, within the host matrix are possible and were considered in our previous work [13].

Fig. 14. (a) Volume fraction dependence of the effective permittivity in the x-direction of the 3D periodic composite displayed in Fig. 11. Inclusions (permittivity "1=30) are ellipsoids of volume fraction f in a matrix of permittivity "2=1. The BIE results are presented for two different values of c: the dashed curve is obtained for a prolate ellipsoid of semiaxes b=c=a/2, whereas the solid curve refers to a prolate ellipsoid of semiaxes b=c=a/4. (b) Same as in Fig. 14 for the effective permittivity in the z-direction.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

411

5.2.1.4. Effect of lattices. Simulations on periodic systems involve two different sets of coordinates: the particle positions in the unit cell and the lattice dimensions. Up to now, we have only considered a simple-cubic (sc) arrangement of the inclusions in the host matrix. Here, we discuss the influence of other types of cubic arrangements: body-centered-cubic (bcc) and face-centered-cubic (fcc). We shall also consider crystalline lattices composed of perfectly conducting spheres (1 ! þ1) of radius a, the permittivity of the host medium being "2 ¼ 1. The volume fraction occupied by the particles can be written as f ¼ n0

4 a3 3 0

ð17Þ

where n0 is the number of particle per unit of internal structure and 0 denotes the volume of the internal structure. Effective permittivity is computed using the BIE

Fig. 15. Notation relating to the unit cell of two-component periodic composite material investigated in the boundary element model computation. The cylindrical inclusion with dielectric constant "1 is periodically arranged in a three-dimensional cubic lattice structure. The dielectric constant of the remaining space is "2 .

412

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

method for the three cubic lattices. The results are tabulated in Tables 1 and 2 together with the solution values by McPhedran and coworkers [21] and Doyle [47]. The first observation that can be made from Table 1 is that the agreement is good. It can be concluded from this agreement that the BIE method has the capability to give high accuracy results.

Fig. 16. (a) Volume fraction dependence of the effective permittivity in the x-direction of the 3D periodic composite displayed in Fig. 15. Inclusions (permittivity "1 ¼ 3) are cylindrical and of volume fraction f in a matrix of permittivity "2 ¼ 1. The full circles are obtained by the BIE method. The solid curve is obtained from Eq. (15). (b) Same for the effective permittivity in the z-direction. Inclusions (permittivity "1 ¼ 3) are cylindrical and of volume fraction f in a matrix of permittivity "2 ¼ 1. The dashed and the solid lines are obtained from Eq. (16) with the apparent permittivities "a="2 and "a=", respectively. (c) Same as in Fig. 16a, when inclusions (permittivity "1 ¼ 30) are cylindrical and of volume fraction f in a matrix of permittivity "2 ¼ 1. (d) Same as in Fig. 16c, when inclusions (permittivity "1 ¼ 30) are cylindrical and of volume fraction f in a matrix of permittivity "2 ¼ 1.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

413

Next we turn to the case of hexagonal lattices composed of perfectly conducting spheres of radius a embedded in a host medium of permittivity "2 ¼ 1. The geometry of the unit cell is displayed in Fig. 17. The number of particles per unitp offfiffiffiffiffiffiffi internal ffi structure is n0=1, and the volume of the internal structure is 0 ¼ ‘2 h 3=2. This type of composite is anisotropic ("x 6¼ "z ) and we observe, from Fig. 18a and b, that the permittivity depends strongly on the specific ratio h‘ for volume fractions f higher than 0.1.

Fig. 16 (continued).

414

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

Table 1 The main characteristics of the different types of cubic lattices sc

bcc

fcc

2

4

Internal structure Cubic array with side 2 1

n0

3

f

a 6

ap

1

fp

  0:523 6

3

2a3 3 pffiffiffi 2 2 pffiffiffi  2  0:740 6

a 3 pffiffiffi 3 2 pffiffiffi  3  0:680 8

The internal structure, the number of particles per unit of internal structure n0, the volume fraction f of the conducting spheres of radius a, the radius ap and the concentration f
corresponding to the maximum packing threshold are function of the type of cubic lattices, i.e. simple-cubic (sc), body-centered-cubic (bcc) and face-centered-cubic (fcc).

Table 2 Comparison of our results with the solutions of McPhedran and coworkers [21] and Doyle [47] for the effective permittivity of different types of cubic arrays of perfectly conducting spheres f

0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70

sc

bcc

fcc

BIE

McPhedran

Doyle

BIE

Doyle

BIE

Doyle

1.000 1.158 1.333 1.531 1.756 2.018 2.333 2.729 3.264 4.083 5.886

1.000 1.158 1.334 1.531 1.756 2.018 2.333 2.729 3.261 4.080 5.887

1.000 1.158 1.334 1.531 1.756 2.018 2.333 2.729 2.363 4.083 5.891

1.000 1.158 1.333 1.529 1.751 2.003 2.292 2.631 3.036 3.531 4.166 5.027 6.344 9.022

1.000 1.158 1.333 1.530 1.751 2.002 2.292 2.631 3.035 3.532 4.166 5.028 6.340 9.023

1.000 1.158 1.333 1.529 1.751 2.002 2.290 2.625 3.023 3.503 4.107 4.889 5.969 7.637 10.920

1.000 1.158 1.333 1.530 1.750 2.001 2.290 2.625 3.023 3.505 4.106 4.888 5.972 7.643 10.925

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

415

5.2.2. Absorbing configurations, scaling and critical exponents Because nonabsorbing spherical inclusions are the exception in nature rather than the rule, we now focus on the case of absorbing particles, with regular geometrical shapes (with boundaries specified by simple functions), dispersed in a nonabsorbing matrix. Such composites may represent a realistic model for many situations in the field of atmospheric science [128]. 5.2.2.1. Spherical inclusions. We begin this presentation by considering a system in which spheres are spatially distributed at the nodes of a simple cubic lattice. In Fig. 19, results are given for the set of permittivity component values: "1=38.2  i22.1 and "2=1.0  i0.0. These values were chosen for the purpose of comparison with the results of Calame et al. [129]. In Fig. 19a and b, we exhibit the real and the imaginary parts of the effective permittivity as a function of the volume fraction f of the inclusion phase. The BIE results compare very well with the data of Calame et al. [129]. There are important things to note here, related to the dramatic alteration of the permittivity for a volume fraction f ffi 0.5 which is close to the concentration corresponding to the touching spheres condition for an infinite simple cubic lattice of uniform size spheres. Another observation one can make from Fig. 19a and b is the following: both "0 and "00 turn down sharply as the volume fraction of inclusions drops below f
and the system shows insulating behavior. We shall defer the interpretation of the conduction threshold volume fraction until the end of this section. We continue by quantitatively analyzing the permittivity data as a function of composition for another set of permittivity component values: "1=80  i106 and "2=2  i0.0, see Fig. 20a and b. The composition dependence of "0 and "00 as f ! f
is presented in Fig. 21a and b. The next question we need to address is whether the transition is governed by universal exponents. In order to understand the divergence

Fig. 17. Hexagonal lattice structure.

416

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

of "0 and "00 , close to f
, in a pragmatic sense we can use a scaling analysis and define (with a slight abuse of notation) two exponents s* and t* as:

s
ð18aÞ "0  f  f
; and

t
"00  f  f


ð18bÞ

where  means ‘‘ asymptotically correct close to f
. The best scaling collapse of the "0 and "00 data determined the values of s* and t*, respectively s*=0.65 0.15 and t*=0.80  0.05. The error in determination of these

Fig. 18. (a) Volume fraction dependence of the effective permittivity in the x-direction of the hexagonal lattice (Fig. 17) composed of perfectly conducting spheres embedded in a host medium of permittivity "2=1. The results obtained by the BIE method are presented for two values of the specific ratio h=‘: solid curve (h=‘ ¼ 1:5) and dashed curve (h=‘ ¼ 2). (b) Same as in Fig. 18b for the effective permittivity in the zdirection.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

417

values is estimated to about 10–15%, again mainly due to the various choices for the volume fraction ranges where data were fitted to a straight line. 5.2.2.2. Cylindrical inclusions. A series of calculations was also performed on rodlike inclusions. More specificaly, the rod radius is kept fixed and the volume fraction of the inclusions changes by varying the rod length l. In our calculations, we varied the rod radius r from 0.05 to 0.45. In Fig. 22a and b the real and imaginary parts of the

Fig. 19. (a) The real part () of the effective permittivity is shown as a function of the volume fraction f of the inclusion phase. Spherical inclusions (permittivity "1=38.2  22.1) are placed in a host matrix material of permittivity "2=1.0  i0.0. Simple cubic lattice. Also shown is the Calame et al. [129] prediction (&). (b) Same for the imaginary part of the effective permittivity.

418

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

relative complex permittivity are plotted versus the volume fraction f, respectively. The most distinguishing feature of the data is the sharp increase in permittivity for f ffi f
, which is particularly evident in the imaginary part. This ‘‘conduction threshold’’ marks the transition from insulating to conducting behavior and is identified with the asymptotic variations of "0 and "00 near f
. We found that the conduction threshold shifts to greater specimen radius as the rod fraction increases, i.e. f
 r2, as can be seen in Fig. 23. Observe that the quadratic law corresponds to the close packing condition since f
occurs when the cylinders touch. It is rather straightforward to quantify this argument. For this system in which one cylinder

Fig. 20. (a) The real part () of the effective permittivity is shown as a function of the volume fraction f of the inclusion phase. Spherical inclusions (permittivity "1=80  i106) are placed in a host matrix material of permittivity "2=2  i0.0. Simple cubic lattice. (b) Same for the imaginary part of the effective permittivity.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

419

Fig. 21. Real part () of the permittivity "0 as a function of the scaling variable f – f
, log–log plot below f
with slope s=0.650.15. The conduction threshold is f
=0.52 0.02. Inclusions (permittivity "1=80  i106) are placed in a host matrix material of permittivity "2=2  i0.0. Simple cubic lattice, (b) same as in Fig. 21a for the imaginary part "00 . Log–log plot above f
with slope 0.80 0.05.

420

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

Fig. 22. (a) The real part of the effective permittivity "0 is shown, in log–log plot, as a function of the volume fraction of the inclusion phase. Rodlike inclusions (permittivity "1=80  i102) are placed in a host matrix material of permittivity "2=2  i0. Simple-cubic lattice configuration of the two-component periodic composite material investigated in this study. Symbols are: (&) r=0.05, (^) r=0.10, (~) r=0.20, (X) r=0.30, (*) r=0.40, (+) r=0.45. the solid lines are guides to the eye. The case of spherical inclusions is also shown (*). The direction of the applied electric field is parallel to the rod axis. (b) Same as in Fig. 22a for the imaginary part of the effective permittivity "00 . Inset indicates the configuration studied.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

421

occupies each unit cell, with L the side length of the unit cell, the volume fraction of 2 inclusion has the form f ¼  Lr 3l. For cylinders arranged on a simple cubic lattice, the
close packing  r2  condition occurs for l=L if the radius is kept fixed. Then f is given
by f ¼  l . Comparing the graphs in Fig. 22a and b, one can observe that the volume fraction dependence of the permittivity is strikingly different for rodlike inclusions from its behavior for spherical inclusions. It is useful to compare the permittivity calculated for a composite based on rodlike and spherical inclusions. Another relevant feature of our simulation results is that "00 is related to f through the heretofore unobserved   scaling behavior "00 ¼ rF ff
: all the data for the imaginary part of the permittivity can be scaled onto a single master curve, see Fig. 24. This scaling demonstrates that, despite the wide variation in the aspect ratio, a single parameter controls the dielectrical behavior of these composites. As of yet, the question concerning the functional form of F(x) is not settled. Of interest is the evolution of the power-law fits applied to the real and imaginary parts of the permittivity in the range 103 < j f  f
j < 101 . Fig. 25a and b show "0 and "00 as a function of the distance from the conduction threshold, ff
, in a double-logarithmic scale. Linear fits to the log-log plots of "0 and "00 give scaling laws which can be characterized by two critical exponents s* and t* according to Eqs. (18a and b), with s*=0.66  0.10 and t*=1.30  0.19 as parameters. Note that these values of s* and t* do not coincidate neither with the universal values of the standard percolation model, nor with the mean field exponents calculated by several authors [36].

Fig. 23. Conduction volume threshold concentration f
as a function of the rod radius, r, of the cylinder.

422

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

5.2.2.3. Ellipsoidal inclusions. We now use simulation data to quantify variations of effective permittivity for ellipsoidal inclusions. To the best of our knowledge, we do not know of any published work on the critical aspects of ellipsoid composite electrical behavior. Fig. 26a and b show the real and imaginary components of the composite material with the long axis of the inclusion parallel to the direction of the applied uniform electric field, and where c is the length of the principal symmetry axis and a is the length of each axis perpendicular to the principal symmetry axis. Here we assume c> a=b which describes a prolate spheroid. Fig. 26a and b reveal the following things. First, at very low loading of inclusions the composite permittivity remains at the level of the permittivity of the unfilled matrix, "2. As the inclusion loading is increased, in a certain region the permittivity increases rapidly by many orders of magnitude. Second, this sharp nonlinear variation in "0 and "00 is sensitive to the aspect ratio (short axis length/long axis length) of the inclusions. However despite their differences all the data do have some features in common. For instance, we have analyzed the scaling of the conduction threshold concentration with the length a. Specifically our data show a square law dependence, i.e. f
 a2, as can be deduced from Fig. 27. One can see that this result is consistent with the close packing condition for a system in which one prolate spheroidal inclusion occupies each unit cell of side length L. The same analysis tells us that the volume fraction of 2 c inclusion is f ¼ 4a . For prolate spheroids arranged on a simple cubic 3L3  2 lattice, the close packing condition occurs when c ¼ L2 , corresponding to f
¼ 6 ac . The data exhibit an unexpected scaling behavior (see Fig. 28) that transcend the large

Fig. 24. The imaginary part of the effective permittivity "00 as a function of the normalized volume concentration f/f
. Same symbols as in Fig. 22a.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

423

Fig. 25. (a) Real part of the permittivity "0 as a function of the scaling variable f – f
, where f f
. In a log–log plot the points fall on a straight line with slope t*=1.30 0.19.

424

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

Fig. 26. (a) Real part of the effective permittivity "0 is shown, in log–log plot, as a function of the volume fraction of the inclusion phase. Ellipsoidal inclusions (permittivity "1=80  i102) are placed in a host matrix material of permittivity "2=2  i0. Simple-cubic lattice configuration. Inset indicates the configuration of the two-component periodic composite material studied. Symbols are: (^) a=b=0.05, (~) a=b=0.10, () a=b=0.20, (*) a=b=0.30, and (*) a=b=0.40. The solid lines are guides to the eye. The case of spherical inclusions is also shown (*). The direction of the applied electric field is parallel to the principal symmetry axis. (b) Same as in Fig. 26a for the imaginary part of the effective permittivity "00 .

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

425

 variations: all the data can be collapsed onto a single master curve, "00 ¼ aF ff
, indicating of a universality in the conductivity property. To analyze the critical behavior more quantitatively, we find that it convenient to use the power-law analysis. The results are shown to scale onto single curves like Fig. 29a and b. As the critical point is approached the numerical data and a least-square fitting provides the exponents s*=0.41 0.07 and t*=0.57  0.04. We observe that they are different from the values obtained for spherical and rodlike inclusions. We will return to a discussion of these values later. From the above, we observe that a parameter that can be manipulated to shift f
is the aspect ratio, i.e. rl, or ac. Note that the square law scaling has been found also in [110] in the context of carbon fibers polymer composites. However, care must be exercised in arriving at any practical conclusions from our results which have been derived for an idealized computing system limited by simple physical considerations, i.e. periodic heterostructures. This result demonstrates the crucial role which the geometry and associated symmetry of the inclusions play in governing the electric field distribution in the system. We have argued in [18,19] that, caution must be exercised if one is to attempt to compare these numerical findings with data obtained in experimental results. For example, it requires a good deal of care to obtain clear experimental results for volume fractions near to the conduction threshold. A major cause of this difficulty is that the conduction transition between the insulating state and the conducting state may be not perfectly sharp, i.e. the rounding for f ffi f


Fig. 27. Conduction volume threshold concentration f
as a function of the length of the semi-minor axis of the ellipsoid a.

426

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

makes it difficult to determine the threshold volume fraction f
accurately. Further investigation of the meaning of scaled data for these systems is warranted. 5.2.3. Multilayered structures As part of our program in characterizing the dielectric properties of heterostructures we have targeted several geometries as structural models for the design of composites with stratified particles based upon the assembly of multilayered structures. In this subsection, we focus our attention on the unit cell which is illustrated in Fig. 30, where the host matrix (permittivity "2 ¼ "02  i"002 ) is filled with inclusions (permittivity "1 ¼ "01  i"001 ) surrounded by an interfacial layer (permittivity "3 ¼ "03  i"003 ). The composite material is a 3D simple cubic lattice arrangement of this unit cell. The uniform external electric field is applied along the principal symmetry axis of the three-phase confocal spheroid-background mixture (a> b=c, with a being the semi-axis of the core inclusion) as displayed in Fig. 30. Within the quasistatic assumption, the largest dimension, a or b, should be smaller than the wavelength of the electromagnetic wave to avoid scattering effects. We denote by f1 and f3

Fig. 28. (a) Real part of the permittivity "0 as a function of the scaling variable f  f
, where f f
. In a log–log plot the points fall on a straight line slope t*=0.570.04.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

427

the filler and the layer volume fractions, respectively. The matrix volume fraction is f2=l(f1+f3). We begin by considering the effect of the permittivity of the interfacial layer on the effective permittivity of the composite material. To illustrate The predictions of the numerical approach proposed above, we take "2=4  i102, "1=20  i2 and "3="03  i0

Fig. 29. The imaginary part of the effective permittivity "00 as a function of the normalized volume concentration f/f
. Same symbols as in Fig. 26a.

428

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

with 1 < "03 < 40, i.e. the interfacial layer is a pure dielectric, as examples. The filler and the layer volume fractions are f1=0.1 and f3=0.04, respectively, corresponding to e/l=0.034, e being the layer thickness of the spheroid. We note here a particular feature for spheroids: because the stratified particle is formed from confocal spheroids, then the interfacial cannot  layer   be uniformly thick. In such a case the quantity e defined as ff31 ¼ 1 þ ae 1 þ be 1 þ ec  1 in our notation represents an average thickness of the interfacial layer. A series of calculations has been made by varying the ratio a/b of the length of the principal symmetry axis to the length of the axis perpendicular to the principal symmetry axis, from 1, corresponding to a sphere, to 1.9, corresponding to a prolate ellipsoid. In Fig. 31a we have plotted the real part of the effective permittivity, "0 as a function of "03 . Two salient features can be observed in this figure. On the one hand, increasing "03 significantly (by a factor of 40) results in a very weak dependence of "0 . On the other hand, changing from a sphere to a prolate ellipsoid has for effect to increase "0 for a given value of "03 . In Fig. 31b, we have compared our results with those derived from Steeman and Maurer’s analytical expression [130] for the same set of permittivity values. Although the data show similar trends, we note that the numerical values differ significantly for spheroids: more explicitely, we find that the predictions using Steeman and Maurer’s approach are higher than our numerical results. We believe that this discrepancy results from a fundamental shortcoming of the Steeman and Maurer equation: it does not directly account for the effect of the multiple interactions between the particles—a consequence of the mean-field treatment of the system. Thus, any equation based on the EMA that excludes multipole interactions is certain to fail over the entire range of volume fraction of inclusions. Next Fig. 32 shows the imaginary part of the effective permittivity, "00 , as a function of "03 . Inspection of the data of Fig. 32 reveals that "00 increases rapidly, then decreases slowly as "03 increases. Another interesting aspect of Fig. 32 is that "00 increases as a/b increases. We take the view that the most

Fig. 30. The unit cell of the three-phase periodic composite material composed of spheroidal (ellipsoids of revolution) inclusions of permittivity "1 and "3, lying in a background of permittivity "2.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

429

distinguishing feature of the data, i.e. "002 < "00 "001 , reflects a screening effect of the interfacial layer. 5.2.4. Comparison with experiment In line of the earlier discussion, experiments have been recently carried out to test the predictions of the numerical model. These experiments have extended the range

Fig. 31. (a) The real part of the effective permittivity "0 is shown as a function of the real part of the effective permittivity of the interfacial layer, "30 . Simple cubic lattice configuration. Symbols are: (&) spherical inclusions; ellipsoidal inclusions (^) a/b=1.1, (~) a/b=1.3, (x) a/b=1.5, (*) a/b=1.7, (+) a/ b=1.9. The solid lines are guides to the eye. (b) The real part of the effective permittivity "0 is shown as a function of the real part of the effective permittivity of the interfacial layer, "30 calculated using Steeman and Maurer’s formula [130]. Same symbols as in Fig. 31a. The solid lines are guides to the eye.

430

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

of available data necessary for testing numerical models [20]. The basic idea of our experiment is to consider a ‘‘pseudoperiodic’’ anisotropic two-component composite material. Here, we investigate model systems in which it is possible to introduce and control the shape anisotropy of inclusions in a host (isotropic) matrix. To provide the simplest geometry for the dielectric characterization of the system, we have manufactured material samples consisting of a cube of dimensions 303030 mm in which N=9 parallel nonoverlapping identical cavities in the form of circular cylinders have been periodically spaced in the plane transverse to the cylinders. The two axially symmetric configurations investigated in the present study are illustrated in Fig. 33a and b. They differ only by the orientation of the cylinders’ axis with respect to the direction of the electric field vector, i.e. respectively parallel or perpendicular. The former sample we call z aligned since the electric field will be applied parallel to the axis of the cylinders, which we define as the z-axis. The latter sample we call x–y aligned. The host matrix is made from poly(vinyl chloride) (PVC). The inclusions contain either de-ionized water, or air. The former (resp. latter) has a higher (resp. lower) dielectric constant than the permittivity of PVC. Varying the volume fraction of inclusion, f, can be realized, either by modifying the radius of the circular cylinder (04r48 mm), or by adding spacers of nominal thickness of 2 mm (see Fig. 33c). The test cell consists of two parallelipipedic copper electrodes with an active section being a square of dimensions 3030 mm. A spring is placed on the upper

Fig. 32. The imaginary part of the effective permittivity "00 is shown as a function of the real part of the effective permittivity of the interfacial layer, "30 calculated using Steeman and Maurer’s formula [130]. Same symbols as in Fig. 31a. The solid lines are guides to the eye.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

431

electrode to maintain a good contact between the electrodes and the sample, in order to avoid the presence of air bubbles manifesting as a parasite capacitance. Capacitance and conductance were measured by means of an impedance analyzer HP 4194A at a frequency of 1 MHz and an amplitude of 3 V/cm, and expressed as "0 and "00 based on the external dimensions of the samples and electrodes [20]. All

Fig. 33. (a) Geometry of the z-aligned two-component composite sample: the electric field is parallel to the principal axis of the non-overlapping circular cylinders, (b) geometry of the x–y aligned two-component composite sample: the electric field is perpendicular to the principal axis of the non-overlapping circular cylinders, (c) same as in Fig. 33a, albeit including spacers in order to vary the volume fraction of inclusions whose radius is kept fixed.

432

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

measurements were performed at room temperature. We have to point out some possible limitations of the experiment. The first is that this experimental model is not rigorously periodic since it is of finite size. Thus it is important that the number of inclusions should be sufficiently large. Only in this case can we assume that this experimental model is quasi-periodic. The second limitation is the actual sensitivity of the measured capacitance to air bubbles due to incomplete filling of the cylindrical inclusions. As the volume concentration of fluid becomes large, this effect is expected to be important. The predictions of the BIE model are now compared with the experimental data. We first carried out some calculations and measurements of the permittivity for z-aligned samples (Fig. 34). A comparison of the calculated real part "0z and imaginary part "00z of the effective permittivity versus volume fraction of inclusions with experimental data, is shown in Fig. 34a and b for inclusions containing voided materials (inclusions with air). It is worth observing that the calculated results can be quite different from the experimental values, for example, Fig. 34b. It can be also seen that the difference between calculated and measured values increases as the volume fraction of inclusions increases. It is observed that "0z and "00z decrease gradually as the void volume fraction increased with numerical values being smaller than the experimental values, over the entire volume fraction range studied. How does shape anisotropy affect the complex effective permittivity? To understand this, we consider first x–y aligned samples. Figs. 35a and b compare the calculated real part "0x and imaginary part "00x of the effective permittivity versus volume fraction of inclusions with experimental data for the voided material. A close look at the values of "0x and "00x shows that there is still some difference between calculated and measured data. Two important features are worth noting. First, we note here that the curves "0x ð f Þ and "00x ð f Þ look very similar to the dependences "0z ð f Þ and "00z ð f Þ. Second, it is instructive to compare these results with those obtained previously. It can be seen that the values of the permittivity along the direction of the electric field are larger than those in the orthogonal direction, and this effect is well outside the deviations of the measurements. This observation is consistent with the theoretical approach given in [130] and the prior experimental work given in [24] on aligned carbon fibers within a polymeric matrix. As the earlier discussion illustrates we have presented a detailed confrontation of experimental and numerical results, in which controlled anisotropy has been produced by selecting a specific geometry with high symmetry. In most cases the qualitative comparison turns out to be reasonably good. Since the correctness and accuracy of the predictions obtained from our algorithm were tested favorably against the results of several analytical and numerical calculations, we believe that the quantitative differences originates from several limitations of the experimental model. Here, we discuss possible explanations for the disagreement between the numerical and experimental data, such as that shown in Figs. 34 and 35. While the cause of this disagreement is not completely understood at this time, the authors can only point out a few reasonable explanations potentially accounting for the sources of experimental error. As mentioned previously, one possible source of disagreement is the necessary finite size of the experimental model. An important question is

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

433

under what conditions is it fair to assume a periodic configuration for the system displayed in Fig. 33. All experimental data presented here were carried out using a reasonable number N of cylinder inclusions but it should be noted that extending the range of the experimental model to a higher number is in principle possible even if the range of N that can be studied is limited by the model. It seems that a low value

Fig. 34. (a) The real part "0z of the effective permittivity as a function of the hole volume fraction, for zaligned samples: numerical simulation (^), experimental (~). Room temperature. (b) Same as in Fig. 34a for the imaginary part "00z of the effective permittivity.

434

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

of N of the experimental model will have an overriding importance at high concentration of inclusions, which is evidenced in our measurements. We have also conducted several experiments with inclusions containing de-ionized water [20]. In many cases the differences are small and can be safely neglected, for example, Fig. 4a of [20]. In some cases we note, however, that the measurements depart significantly from expectations based on our numerical model, for example, Fig. 6a of [20]. A

Fig. 35. (a) The real part "0x of the effective permittivity as a function of the hole volume fraction, for x–y aligned samples: numerical simulation (^), experimental (~). Room temperature. (b) Same as in Fig. 35a for the imaginary part "00x of the effective permittivity.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

435

possible source of disagreement between the simulation and experimental data may be due to the presence of air bubbles caused primarily by a less than satisfactory filling of the inclusions. It seems likely that this could explain why the nominal values of "0x and "00x are lower than those of "0z and "00z (compare Figs. 34 and 35). Moreover, and interestingly enough, we want to point out that ill-controlled filling of the inclusions does have a quantitative impact on permittivity which is consistent with the lower (resp. higher) observed "0 and "00 than the corresponding calculated values of "0 and "00 for de-ionized water (resp. air), for both x–y and z aligned samples. The former case corresponds to a material medium filling the inclusion having higher "0 and "00 than the "0 and "00 of the host matrix while the latter case corresponds to the opposite behavior. 5.3. Discussion and final remarks Some comments seem in order at this point. The ab initio calculations reported here provide high precision data for analyzing the permittivity of dielectric heterostructures and highlight a number of complexities not previously noted for these systems. They also make several predictions that may be tested by the experiment which may be very useful for confirming or discarding a theoretical study. Our approach to modeling the effective permittivity is largely based on the QA that is expected to be accurate at relatively low frequencies. The outcome of the simulations carried out in the present work is a quantitative measure of the effects of changing the volume fraction, the shape, and the orientation of the inclusions within the host matrix on the effective permittivity of two-component heterostructures. Several 2D and 3D configurations of experimental interest have been treated. The derivation presented in this paper has also the advantage of showing clearly the underlying unity of the subject matter.We can draw a few immediate conclusions from these results. Overall, calculations are very good, particularly when it is taken into account that the computational cost expended in obtaining them on a relatively small computer and the memory requirements are low. We believe that by now we have accumulated extensive numerical evidence to support our view that our BIE approach is exact for periodic arrays of inclusions embedded in a matrix. The remarkable agreement of our data with analytical calculations and experiment increases the confidence that can be placed in that numerical procedure. We now return to the issue of power law dependence for the effective permittivity in the periodic models considered here and the question of universality. Collectively, the data reported here demonstrate that the shape anisotropy of inclusions in dielectric heterostructures affect dramatically the effective complex permittivity. Remarkably, we found that the square law dependence of the conduction threshold versus the aspect ratio applies for both rodlike and ellipsoidal inclusions, over a wide range of aspect ratios. Based on the preceding observations it is important to make a number of qualifying remarks. First, the effective permittivity of dielectric heterostructures composed of inclusions having shape anisotropy can be an extremely sharp function of volume fraction of inclusions, depending on the aspect ratio, and the orientation, and the conductivity of the inclusions. The results we have

436

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

obtained are consistent with previous observations made by several groups [108]. Second, one of the most interesting conclusions which arises from the above results is that the transition from insulating to conducting behavior is critical, exhibiting scaling laws. There are indications of unconventional critical exponents from powerlaw volume fraction dependences of the permittivity. Study of the scaling laws reveal lessons for how to make dielectric heterostructures having a low conduction threshold. A question that may arise upon evaluating these data concerns the actual nature of this transition: if scaling and universality exist in the critical region, to which universality class does the transition belong? It would be highly desirable to have this question answered on solid theoretical ground. There is considerable value in studying the features of such simple systems to determine what, if any universal characteristic emerge, and how these can be manipulated through the geometrical shape of the inclusions. Third, our results illustrate the importance of properly describing the alignment and contact of inclusions in the matrix, a result that opens up interesting possibilities for the future analysis of experimental data. Shape anisotropy and spatial orientation in the dielectric properties of heterostructures are of considerable interest, especially in the light of the suggestion that the conduction threshold is related to the actual area of contacts between the inclusions. The essential point is that the different behavior between isotropic and anisotropic inclusions must reflects the surface area of the inclusions. What general statements can we now make about the connection between conduction threshold and contact? The answer can be guessed on symmetry grounds. To explain, we may state matters as follows: Spherical shape provides a minimum of surface area for a given volume fraction and a fixed number of inclusions. Similarly elongated shapes, which deviate from spherical, provide less surface area for a given number and volume fraction of inclusions. Consequently, spherical inclusions must occupy a higher volume fraction before guaranteeing contact for the existence of a continuous path that crosses the sample. We may observe that, for inclusions having anisotropic shapes, the full rotational symmetry of the problem is broken at the critical point, a unique preferred direction is picked with different critical exponents than the isotropic problem: this is apparent from a comparison of s* and t* in Table 3. As this work indicates, the sharp increase in the permittivity when plotted as a function of volume filling fraction characterizes higher electric multipole interactions because all of the inclusions contact at once. We note also that the values of f
are consistent with analytical results from theoretical predictions based on the concept of excluded volume associated with 3D objects [108]. The excluded volume is much larger than the actual volume for extended shapes, such as rods and ellipsoids. Then, the resulting f
may be lower than for an array of spheres. For spheres the actual and excluded volumes are the same. Equally significant, we observed a strong orientational effect on the permittivity [19]. This again indicate that the more extended shape which deviates from spherical provides what is needed for contact: they are characterized by more favorable depolarization factors. It is equally true, however, that physical mechanisms other than contact are also involved if we are not dealing with an idealized computer model but with real engineering materials. Many questions still remain as to what factors dominate, particularly when the interfacial

437

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

contact is mediated by adsorption zones that introduce parameters such as adhesive interactions at contact, tunneling phenomena, and partial wetting. This lack of knowledge about structures has limited an understanding of the physical properties of these materials. For random systems, we expect that two factors govern the contact between inclusions: geometrical and statistical. Geometrically, the number of contact points per inclusion is related to the surface area of the inclusion; statistically the contacts involving a given component increase with the increase in the amount of that component. Even though an acceptable description of the complex 3D structure with the interface surfaces oriented in random directions with respect to the applied electric field has been reported, it does not provide a satisfactory description of the dielectric properties of real heterostructures for which specific interactions between the inclusions and the host matrix need be taken into account. Moreover, in (random) composites containing complex-shaped constituents, an exact calculation of higher multipole interactions is intractable. Turning finally to the experimental verification of these numerical calculations. The results obtained by this method were tested experimentally on a specific geometry of inclusion involving two systems (inclusions filled either with air, or deionized water within a polymeric matrix) [20]. While the overall features in the experimental results agree with the simulation, there were still some discrepancy between them. There are at least two important factors that contribute to these differences: (1) finite-size effects in our experimental results, and (2) incomplete filling of the inclusions. These two features were discussed in detail in [20]. Additional experiments would appear to be required to fully understand the polarization mechanisms that control the permittivity of heterostructures.

Table 3 Results of fitting the real and imaginary parts of the effective permittivity of two-component composite materials to a simple power law giving the critical exponents s* and t* Inclusion

s*

t*

Sphere Rod Ellipsoid

0.650.15 0.660.10 0.410.07

0.800.05 1.300.19 0.570.04

Inclusions of permittivity "1=80  i102 are placed in a host matrix of permittivity "2=2  i0. Simple-cubic lattice configuration. For a closely packed set of spheres randomly arranged in a matrix, Batchelor and O’Brien [Proc R Soc Lond A 1977;355:313] used lubrication theory and obtained for perfectly conducting spheres a logarithmic divergence of the effective thermal (or electrical) conductivity of the material, i.e. in our notation s*=0. However, this result is by no means conclusive since it takes into account only dipolar interactions. See also [36] for a more recent approach. Here, the author has derived analytical expressions for the divergence of "0 and "00 of cylindrical inclusions at the critical point. In this approach, the critical exponents s*=1 and t*=1 are mean-field exponents. However, as has been widely discussed in the literature, the predictions of mean-field theories often contradict many of the experimental data on random composite materials because these theories do not describe accurately the actual multipole interactions at the critical point.

438

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

6. Conclusions and future directions The basic goal of the research described in this article is an effort spent on a better understanding of the dielectric properties of heterogeneous materials via computational electromagnetics. A real experiment on real composite materials is a very powerful tool. However, the results are necessarily an average over mesoscopic to macroscopic length scales. As has been realized independently many times, the complexity of the many-bodied interactions and the heterogeneous environment experienced by the individual particles in matrix-particles composites is generally not amenable to a simple analytical description. Much more details can be probed in a computational experiment. Computer simulation aids in understanding the underlying fundamental physics of these mixed media as well as is useful in engineering applications. In this article we have given the reader a mere glimpse of what has been done with a fully ab initio numerical algorithm that is capable of producing exact calculations of the dielectric characteristics in model (translationally invariant) materials. Let us note once more that the underlying quasistatic assumption inherent in our derivation necessarily restricts the range of particle size and typical separation between inclusions that can be studied. Continuum electrostatics provides a natural and intuitive description of a great many of real composite materials, because their granularity is not resolved by commonly achievable wavelengths. A number of interesting results emerge from the foregoing calculations which can be subject to experimental test. Taking all the current results together, it can safely be concluded that these lattice simulations indicate that we are able to capture the main aspects of the dielectric behavior consistent with the assumptions made in the previous sections. As a motivation mentioned in the Introduction, this understanding may lead to better models of the real experiments, suggest methods of data analysis, be used to benchmark and validate computer models, etc. Thus numerical studies of model materials complement experimental studies of real materials. In many instances, the dielectric behavior is controlled by the geometry and size of the separate inclusions in the matrix. Through extensive benchmarking by using periodic geometries, it has been found that the BIE approach can give very good results as far as we are concerned with very long-wavelength physics. Hence the basic definitions and the the numerical procedure were discussed. We also showed that the computed properties can be compared fairly well with their experimental counterparts for a dielectric composite model. The purpose of our work was not to decide which EMT are good or bad, in general, but to find which of them can be expected to lead to reasonable results when applied to the considered case of composite material. We hope that our discussion has demonstrated that our treatment of periodic composite materials is quite reliable. Overall, the results suggest that it will be always difficult to reproduce the experimentally measured dielectric characteristics, when limited to dipolar terms and the problem would be even harder if one wanted an approximation which worked for currently available (random) composites. We have only been able to touch on some of the unique dielectric properties of heterostructures. The BIE algorithm is sufficiently general so that heterostructures with a wide range of topological properties can be modeled through proper specification of parameters such

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

439

as permittivity of the constituent materials, size, geometrical shape, and topological arrangement in the heterostructure. We invite the reader to look at the published articles of our group cited in the references for other results [13–20]. As we commented earlier, the technique that we developed can also be used to describe other physical properties of ordered composites, for example, the equations governing the electromagnetic behavior of these systems are formally equivalent to the heat transport problem [131]. Prediction of the effective permittivity of (real) composite materials from first principles is a largely unmet goal of theoretical solid-state physics though a few groups have reported such calculations. We have learned through empirical studies, analytic and computational approaches that our understanding of the polarization mechanisms in complex condensed matter systems is still far from complete and is capable of surprises. But the body of computational evidence that we have outlined here fits well into the actual trend towards computational physics approaches of heterostructures. Yet we believe that these results are only the tip of the iceberg. We have uncovered some generic computational principles, but the answers to the open questions could hide new concepts and ideas that might turn out to be just as exciting as those we have encountered so far. Many important fundamental questions remain unanswered such as: To which universality class does belong the critical behavior of the conduction threshold? Is it related to the spatial dimension, or to the type of heterostructure? In order to understand the physics beyond the QA, one needs a better understanding of the issue of grain-scale description of strong multiple scattering. When do inclusion boundaries become significant, and what is their effect? When particles are not precisely identical, but have a range of radii (and possibly shapes, etc), what is the effect on the effective permittivity? All of these exciting questions are topics of current research and are apparently far from being settled. Addressing these and similar questions is important for obtaining clues about the dielectric properties of composite materials. It is our intention to report on these computational aspects in a forthcoming publication. Understanding these things will open the way to designing new engineering materials where the dielectric properties can be tailored to a particular application with a high degree of precision At this point it seems appropriate to consider what directions research in this field will take in the next few years [132]. In our opinion, understanding the physics of the charge and electromagnetic wave transport in complex random heterostructures is one of the key problems in condensed matter. This problem has so far been understood only in the mean-field level [3,33,48]. However, EMT have the disadvantages typical to all mean-field theories since it ignores the fluctuations in the system. It assumes that the local electric and magnetic fields are the same in the volume occupied by each component of the composite material. In some cases the local field fluctuations cannot be ignored, i.e. near the percolation threshold. Composite materials are prototypical of correlated systems where, many-body effects, resonance effects, and interface effects can be at play simultaneously, and where usual simplifications that neglect some interactions to study others in detail can not work. We are convinced that the comprehension of complex problems posed by systems which exhibit the effects of inhomogeneities has the potential to produce a

440

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

substantial advance in many areas of condensed matter physics, and will eventually shed light into other unsolved problems, for example, to name for a few, inhomogeneities arising from crystal imperfections, bonding defects, or impurity clusters in semiconductors, inhomogeneities in materials which undergo metal-semiconductor transitions, and inhomogeneities in disordered systems such as glasses. From our personal perspective, we clearly see a multitude of scientific challenges. The immediate focus is likely the finite-size numerical study of heterostructures with broken-translational symmetry. The above described approach contains a great deal of the dielectrics of heterostructures without the complication of the irregular or random geometry. In sum, our numerical findings are interesting, but preliminary; a continuing systematic study of wave transport properties in more general (random) inhomogeneous systems is necessary to develop new principles for creating structures of controlled length scales that would result in novel materials with distinct (and unusual) properties. The simulation of physical systems with a wide range of temporal and spatial scales poses a considerable challenge to the computational scientist and is still outside the realm of physics-based multiscale first-principles approach. This will require going back to the fundamental question of what the percolation transition is. A crucial component to this computational approach involves determining how characteristic physical quantities and lengths, e.g. the correlation length
, describing the dielectric system scale near the singularity. Describing the current state of the statistical description of random heterostructures in great detail is beyond the scope of this paper, but we believe it is reasonable to say that a computational model faithfully representing the entire gamut of the physics underlying the complex dielectric behavior in our systems does not exist at present, and its development is highly desirable. The microscopic local fields near the interfaces may differ considerably from those in the bulk regions, due to the existence of the scattering fields from boundaries. Work is in progress to extend our method to this case, i.e. to describe the interface between the different materials since it plays an essential role in any device action. It is clear that there are many possibilities in terms of models and parameters, and the detailed predictions for the electromagnetic properties will vary according these choices. The level of detail incorporated in these modeling strategies ranges from electric dipole approximation to higher multipole approximation. The successful description of these systems in terms of Maxwell’s equations relies on the correct (long-range) multipolar forms for the field vectors D and E in these equations. Within this context, the importance of largescale parallelization of efficient and space-saving algorithms for the computation of a full microstructural morphology is by now widely recognized. If physically realistic backgrounds could be derived this way, then this might be a satisfactory outcome. Physics is a paradigm where theory and experiment must work together. As the earlier discussion illustrates, there are not many experimental data for comparison with the detailed analytical and numerical data. The lack of an analytical solution to 2D and 3D heterostructures has not curtailed experimental work in the field. Experiments should tell us a lot. Generally speaking, it is recognized that proper quantification of microstructural (or mesostructural) behavior is necessary for the optimization of materials properties. Experimentalists are beautifully tooled up, as

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

441

never before with sensitive local. Most experiments aimed at observing heterogeneity are based upon measuring the response in a spectrally or spatially selective fashion. Depending on the length scale of the material medium under consideration, various experimental techniques, including both electromagnetic wave (m to mm) and electron (nm to m) microscopy have been employed in structural interrogation. Over the past few years several authors have emphasized the wealth of information available from local probes of the conductivity of composite materials. Improvement of AFM has allowed this vision to become a reality. Through direct observation of the local resistance in ultramicrotomed samples of carbon black-filled polymers, a wealth of information can be obtained on the carbon black aggregates and agglomerates distributions. The capability of AFM to achieve high spatial resolution down to the nanometer level, real space imaging of polymer surfaces is now well documented in the literature [133]. This gives images from which microand meso-scale details can be deduced. There are certain problems with the interpretation of the images, although these could be due to effects related to sample preparation. Structural information of this type is of utmost importance to get an unambiguous signature in the data that can be used to make a valid claim for a percolation transition and how it is affected by the processing conditions [63]. Twocomponent microstructure can be represented by field variables, for example, composition and long-range order parameters that distinguish the parent and product phases. Another question appears to be closely connected to the above issue: namely, how one might characterize the randomness from measurements (or computations) of averages and higher moments of composition fluctuations. To elucidate the physics behind these correlations the technique introduced recently by Torquato [6,48] appears as promising. For such systems, progress has been hampered by the lack of detailed experimental results for the interfacial structure. Experimental methods which probe interfaces must inevitably extract surface data admist that from the response of an overwhelming number of bulk particles. The keys to the future are our ability to improve the properties of materials by structuring and coupling them on different length scales and our ability to produce these materials in commercially viable quantities. Research to date has only scratched the surface of the multiscale analysis of composite materials. By this we mean that one set of calculations at a given level is used to evaluate quantities for use in a more approximate or phenomenological computational methodology at longer length scales. It is also important to establish whether the reconstruction method of random media, given limited microstructural information, yields a unique reconstructed solution, i.e. the so-called degree of achieved topological uniqueness [6,48]. Adding another layer of complexity to the problem of the dielectric behavior is the appearance of resonances, for example, Mie resonances, for the different components of the composite. The difficulties plaguing the task of probing the basic mechanisms governing the polarization behavior stem from the complexities of the microstructure that dominates the response of the material. It is for experiment to show theory the ground it must walk on. There is an a lot to learn: our understanding of real composite materials is little more than phyllotaxic and being seriously outstripped by the discovery of even more exotic and artificial ones. While reporting error bars is an

442

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

important part of any experimentalist’s work, theorists should also devote significantly more effort to quantifying the uncertainties of their calculations. Another major theme to emerge from the discussions of future directions is the broad scope for work on problems in biological physics and other chemical systems of practical interest. Specifically, the cross-coupling in interfaces of force fields arising from electrical, mechanical, chemical and entropic potential gradients has captured the imagination of the physics community, exhibiting new features [134]. As a consequence, reliable theoretical models are needed in order to interpret experimental data and also for the diagnostic and therapeutic medical applications [135]. The multiple relaxation processes and the overall complexity of the interfacial processes cause serious problems to device designers when trying to simulate sensitive sensors for biocompatible applications, for example, tissue-based biosensors. Notwithstanding the large body of both experimental and theoretical work the physics of these effects is not yet fully understood and will be investigated in more detail in future publications. So far, our discussion was restricted to effective linear permittivities. The significance of linear response is that the permittivity is independent of the applied electric field. In the context of nonlinear effective dielectric behavior, some progress has been made during the last two decades [3,36,136,137]. A proper theoretical study of the effects of nonlinear dielectric properties of metal/dielectric composites was undertaken by several investigators including Stroud et al. [136], Agarwal and Gupta [137], and Ponte Castaneda [138]. Their theories which have been the basis of most discussions concerning the subject suggest that important gaps remain in our understanding of the properties of this entire class of nonlinear dielectric materials. Future work will be attempted in our group shortly to study the optical response of composite nonlinear materials. As an additional note, we point out that recent advances in fabrication techniques are providing the opportunity to engineer the properties of composite materials by varying the underlying structure on several length scales including the mesoscopic scale. For example, there has been a growing interest in the study of photonic crystals (artificial periodic dielectric composites, with lattice spacing comparable to the wavelength of the electromagnetic wave) which is inspired to a large extent by the pioneering suggestion of Yablonovitch [139] and John [140] that a photonic band gap could lead to inhibited spontaneous emission of molecules and electromagnetic waves localization. The intriguing electronic and optical properties of photonic crystals have unleashed a flood of interest. For instance, it has been shown that passive elements such as filters, optical switches, cavities, and high-efficient light emitting diodes can be substantially improved if constructed on the basis of photonic crystals [141]. The present method of calculation of the effective complex permittivity can be used to reproduce a photonic structure although this is not its intended application. Another line of research involves the physics of wave transport in nanophases. Since Eric Drexler brought the term ‘‘nanotechnology’’ into vogue in his 1986 book [142], reports on nanostructured-materials and nanoscale devices and their unique physical properties which may be useful in various technological applications have

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

443

permeated the literature in the last decade. The retentless drive toward solid state structures of ever decreasing size scales has been a major driving force in microelectronics with research being carried out all over the world. The ultimate goal (‘‘holy grail’’) of nanotechnology is to synthesize devices as small as a molecule or a cluster of atoms. This point was made by Richard Feynman, at the 1959 annual meeting of the American Physical Society, who pounded on the idea that the entire 24 volumes of the Encyclopedia Brittanica could be written on the head of a pin [142]. The smaller the scale, the better the prospects for emulation. There are several encouraging aspects to this pursuit. Small particles in the size range from a few to several nanometers often exhibit structures and properties which are significantly different from those of the corresponding bulk materials. Nanophases are a unique class of materials that interpolate between two extremes, atomic and bulk phases. Interest in nanophases has expanded as investigators have recognized that many of the properties of finely divided matter strongly depend on the interfacial properties of the constituents by virtue of the high fraction of the overall material at or within a few atomic layers of the surfaces or in the vicinity of an interface as well as of the confinement of electrons, excitons, and photons in small volumes. The unusual electronic properties, if exploited, could lead to nanometer-scale electronic devices. Such nanostructures can be expected to have important applications in, for example, sensors, and storing information systems. Despite many advances and massive investment programs, the nature of many unique behaviors remains controversial and is still one of the key issues in current debates in nanophysics. Much more can be done and ultimately harness the electromagnetic properties of nanostructured materials. New research directions are emerging and are being developed aggressively for a variety of applications as well as for accessing new regimes of basic experimental research. One recently recognized class of finely divided matter, that is now in the embryonic stage, concerns magnetoelectric nanocomposites [143,144]. Future studies are aimed at testing how computational electromagnetics approaches should also serve to deepen our understanding of the role of grain boundaries and/ or interfaces to the determination of the dielectric properties of nanostructures. It is expected that this new field will spur the development of new devices which cannot be realized with existing semiconductor-based electronics [145]. These hybrid nanostructures form a fascinating melting pot for studying the interplay between dielectric and magnetic properties, often revealing new and unexpected physics. In fact, the addition of the spin degree of freedom to conventional electronics knowhow is paving the way to the appealing field dubbed ‘‘spintronics’’. This emerging area includes the active control of carrier spin dynamics and transport in electronics materials, with the final goal of nanolithic integration of the nanoscale with larger structures. For example, colloidal crystals, formed as tiny particles which are deposited from suspension in a fluid onto the surface of an electrode, can be controlled by adjusting the strength and frequency of the applied electric field [146]. The widespread natural occurence of multiscale systems provides further motivation for research in this direction for the coming decades. In addition, we suggest that the BIE approach can be used to investigate electromagnetic properties of materials with simultaneously negative permittivity and

444

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

permeability, or left-handed materials (LHM). In 1968 the Russian physicist Victor Veselago [147] theoretically investigated the electromagnetic consequences of a LHM and predicted that such a medium would have dramatically different propagation characteristics stemming from the sign change of the group velocity, including reversal of the Doppler shift and Cherenkov radiation, anomalous refraction, and even reversal of radiation pressure to radiation tension. However, at that time, these effects could not be experimentally verified since, as Veselago pointed out, materials with negative permeability were not available. Recently, Pendry and coworkers and Smith et al. discovered a microwave-plasma thin-wire structure exhibiting negative permittivity below the electronic plasma frequency [148]. LHMs can potentially be applied to novel applications such as antenna radoms, band-pass filters, microwave couplers, and high-resolution optical lithography in the semiconductor industry. Computational electromagnetics techniques have given a new methodology for answering old questions, i.e. the direct and inverse problems that we defined in the Background section. As these and other questions get robust answers, computational electromagnetics will become a toolbox with which to construct new technological devices. This exciting area of research requires expertise from such disciplines as materials science, computattional engineering, condensed matter physics, and chemistry. The trend towards an integrative approach in many physics departments over the past decade has encouraged a multidisciplinary approach. There is good reason to imagine that one could achieve this goal by building theoretical models running computational simulations. But, we must never forget that is there is subtle but fundamental difference between simulating reality and building it. We hope that this work inspires experimenters and theorists to attack some of the many problems in this exciting and challenging field. In short, there are rich opportunities for a continued impact of materials simulation, and the technological capabilities appear to be up the challenge. The future success of any rational design of dielectric heterostructures will depend on the generalizability of fundamental studies such as the present one.

Acknowledgements Our work has been stimulated by discussions with many colleagues and workers in the field. Throughout this paper, we have concentrated on our own work towards dielectric heterostructures, especially on the late computational parts of it. We must apologize to the many whose work we were not able to treat in the depth it deserved in this review, a sin we have tried to atone for by including an extensive bibliography. The work stemming from a research group cannot be better than the persons carrying it out. It is a pleasure to thank warmly B. Sare´ni, L. Kra¨henbu¨hl and A. Boudida, with whose collaboration some of the above work was done in the initial stage of this project. Laboratoire d’E´lectronique et Syste`mes de Te´le´communications is Unite´ Mixte de Recherche CNRS 6165 and Centre de Ge´nie Electrique de Lyon is Unite´ Mixte de Recherche CNRS 5005.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

445

Appendix A In this appendix, we list some interpolation functions associated to triangular boundary elements commonly used in the BIE method for the computational implementation. Assuming that a surface is discretized with triangular boundary elements as described in Fig. 36a–c, interpolation functions are written in terms of triangular area coordinates 1 , 2 . For linear elements: 8 < l1 ¼ 1  1  2 ðA1Þ l ¼ 1 : 2 l3 ¼ 2 where 1 and 2 are defined in Fig. 37. Then, the values of the cartesian coordinates for any point inside the triangular element can be deduced from the interpolating functions and the coordinates of each node: 8 < x ¼ l1 x1 þ l2 x2 þ l3 x3 ðA2Þ y ¼ l1 y1 þ l2 y2 þ l3 y3 : z ¼ l1 z 1 þ l2 z 2 þ l3 z 3 similarly, the potential and its normal derivative are obtained from Eq. (A3). 8 l2 V2 þl3 V3 < V ¼ l1 V1 þ    @V @V @V @V ðA3Þ ¼ l þl þl 1 2 3 : @n @n 1 @n 2 @n 3 For quadratic elements: 8 l1 ¼ 3 ð23  1Þ > > > > l2 ¼ 1 ð21  1Þ > > < l3 ¼ 2 ð22  1Þ l4 ¼ 43 1 > > > > l ¼ 41 2 > > : 5 l6 ¼ 42 3

ðA4Þ

where 3 ¼ 1  1  2 . Interpolation functions for triangular elements of higher order or for quadrilateral elements of any order are the same used in the finite element method and can easily be found elsewhere [123,124].

Appendix B The purpose of this appendix is to present the details of the calculation for the effective permittivity which leads to the result of Eq. (9). Briefly put, we consider an

446

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

inhomogeneous mixture consisting of inclusions of constituent 1 ("1 ; 1 ¼ 1) surrounded by a background material of constituent 2 ("2 ; 2 ¼ 1), submitted to an average macroscopic field E. The average displacement vector D is linearly related to the average electric field. Here we restrict our consideration to cases where "1 et "2 are real and positive, i.e. cases where there is no absorption and free propagation

Fig. 36. Triangular boundary elements: (a) linear, (b) quadratic, and (c) cubic.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

447

Fig. 37. Triangular area coordinates.

exists for each of the two components. The coefficient of proportionality defines the effective permittivity " which can be written as " ¼ "2 þ n

El ; E

ðB1Þ

where n is the number of inclusions per unit volume,  is the polarizability of each inclusion and El stands for the local electric field. To compute the effective permittivity we can apply Yaghjian’s formula [149] of the local electric field:

El ¼ E þ

1 L:nEl ; "2

ðB2Þ

where L is a depolarization dyadic or a source dyadic that depends on the shape of the inclusions. Eq. (B2) is in fact quite general, holding for arbitrary shape of inclusions and their topological arrangement. The obvious question, of course, is what value of a to insert therein. The polarizability  can be deduced from de dipole moment p of each particle: p ¼ El :

ðB3Þ

The evaluation of the dipole moment requires the knowledge of the internal field of an inclusion, i.e. to find the solution of the Laplace’s equation [46,150]. Both, Eqs. (B1) and (B2) are based on the premise that the binary mixture contains no voids. It may be objected that Eqs. (B1) and (B2) cannot be applied directly to a monodisperse mixture of spherical elements, since such mixture cannot fill space entirely, even at closest packing. Bruggeman [82] was well aware of this argument and in order to fill space with spherical elements, he used a heterodisperse set of ever smaller spheres.

448

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

References [1] Landauer R. In: Garland JC, Tanner DB, editors. Electrical transport and optical properties of inhomogeneous media. AIP Conf. Proc. No. 40. New York: AIP; 1978. p. 2 See also Landauer R. J Appl Phys 1952;23:779. [2] Priou A, editor. Dielectric properties of heterogeneous materials. Progress in Electromagnetics Research. New York: Elsevier; 1992. [3] Bergman DJ, Stroud D. In: Ehrenreich H, Turnbull D, editors. Solid state physics, advances in research and applications, vol. 46. San Diego: Academic; 1992. p. 147. [4] Nan CW. Prog Mater Sci 1993;37:1. [5] See the series of papers: Hori M. J Math Phys 1973;14:514, ibidem 14, 1942, 1973; Hori M, Yonezawa F. J Math Phys 2177;15:1974, ibidem 16, 365, 1975; Hori M. J Math Phys 1722;16:1975, ibidem 18, 487, 1977. [6] Torquato S. Appl Mech Rev 1991;44:37; See also Cule D, Torquato S. J Appl Phys 3428;86:1999; Torquato S, Lado F. Phys Rev B 6428;33:1986; Yeong CLY, Torquato S. Phys Rev E 1998;57:495; Sheenan N, Torquato S. J Appl Phys 2001;89:53. [7] Hale DK. J Mater Sci 2105;11:1974. [8] Newnham RE. Annu Rev Mater Sci 1986;16:47. [9] There are a large number of experimental studies of the physical properties of particles filled polymer systems. Progress in understanding the stucture–property relation of filled polymers is summarized in recent reviews such as Donnet JB, Bansal RC, Wang MJ. Carbon black, science and technology. 2nd ed. New York: Dekker; 1993; See also Medalia AI. Rubber Chem Technol 1978;51:437; Wang MJ. Rubber Chem Technol 1978;71:520. [10] Jonscher AK. Nature (London) 1977;267:673; See also Jonscher AK. Dielectric relaxation in solids. London: Chelsea Dielectric Press; 1987; Jonscher AK. Universal relaxation law. London: Chelsea Dielectric Press; 1996; and more recently Jonscher AK. IEEE Trans Elec Insul 1992;27:407; Jonscher AK. J Phys D: Appl Phys 1999;32:R57. [11] Dyre JC. Phys Rev B 1993;48:12511; See also Dyre JC, Schroder TB. Rev Mod Phys 2000;72:873. [12] Erman B, Mark JE. Structures and properties of rubberlike networks. New York: Oxford University Press; 1997; See also Treloar LRG. The physics of rubber elasticity. 3rd ed. Oxford: Clarendon; 1975; Steeman PM. PhD thesis, Technische Universiteit Delft, Delft, The Netherlands, 1992. For an analytic approach to the interfacial polarization in heterogeneous systems; we also refer the interested reader to Fu L, Macedo PB, Resca L. Phys Rev B 1993;47:13818; Hanai T. Kolloid, Z 1960;171:23. [13] Sareni B, Krahenbu¨hl L, Beroual A, Brosseau C. J Appl Phys 1688;80:1996. [14] Sareni B, Krahenbu¨hl L, Beroual A, Brosseau C. J Appl Phys 4560;80:1996. [15] Sareni B, Krahenbu¨hl L, Beroual A, Brosseau C. J Appl Phys 2375;81:1997. [16] Boudida A, Beroual A, Brosseau C. J Appl Phys 1998;83:425. [17] Beroual A, Brosseau C, Boudida A. J Phys D: Appl Phys 1969;33:2000. [18] Brosseau C, Beroual A. J Phys D: Appl Phys 2001;34:1. [19] Brosseau C, Beroual A, Boudida A. J Appl Phys 7278;88:2000. [20] Beroual A, Brosseau C. IEEE Trans Dielectrics EI 2001;8:921. [21] McPhedran RC, McKenzie DR. Proc R Soc London Ser A 1978;359:45; McKenzie DR, McPhedran RC, Derrick GH. ibid 1978;362:211; Nicorovici NA, McPhedran RC. Phys Rev E 1945;54:1996; McPhedran RC, Nicorovici NA. Physica A 1997;290:173.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

449

[22] Shen LC, Liu C, Korringa J, Dunn KJ. J Appl Phys 7071;67:1990; See also Liu C, Shen LC. Model Simul Mater Sci Engin 1993;1:723. [23] Ghosh PK, Azimi ME. IEEE Trans Dielec Insul 1994;1:975. [24] Cukier RI, Sheu SY, Tobochik J. Phys Rev B 5342;42:1990; See also Sheu SY, Kumar S, Cukier RI. Phys Rev B 1431;42:1990; We also refer the interested reader to Lam J. J Appl Phys 4230;60:1986. This article considers the calculation of the effective permeability of a single cubic lattice of identical conducting magnetic spheres as a series expansion in ascending powers of the inclusion volume fraction. [25] Schwartz LM, Banavar JR. Physica A (Amsterdam) 1989;157:230; See also Schwartz LM, Banavar JR. Phys Rev B 1989;39:11965. [26] Sheng P. In: Ericksen JL, Kindeklehrer D, Kohn R, Kons JL, editors. Homogenization and effective moduli of materials and media. New York: Springer Verlag; 1986 See also Ma H, Zhang B, Tam WY, Sheng P. Phys Rev B 2000;6:962. [27] Pekonen O, Ka¨rkka¨inen KK, Sihvola AH, Nikoskinen KI. J Electromagn Waves Applicat 1999; 13:67. [28] Ka¨rkka¨inen KK, Sihvola AH, Nikoskinen KI. IEEE Trans Geosci Remote Sensing 1303;38:2000. [29] Sihvola A. Electromagnetic mixing formulas and applications. London: IEE Publishing; 1999; See also Sihvola AH. IEEE Trans Geosci Remote Sensing 1989;27:403; Sihvola AH, Lindell IV. J Electromag Waves Appl 1989;3:37. [30] Ishimaru A. Wave propagation and scattering in random media. New York: Oxford University Press; 1997; See also Tsang L, Kong JA. J Appl Phys 3465;51:1980; Tateiba M. Radio Sci 1987;22:881. [31] Kunz KS, Luebbers RJ. The FDTD method for electromagnetics. Boca Raton: CRC; 1993. [32] The art of numerical simulation for solving Maxwell’s equations in the time domain has made significant progress over the last decade and a number of variations of the FDTD are to be found in the literature, for example, see Taflove A. Computational electrodynamics—the finite-difference time-domain method. Norwood, MA: Artech House Inc; 1995. [33] Ericksen JL, et al., editors. Homogenization and effective moduli of materials and media. New York: Springer; 1986; See also Dal Maso G, Dell’Antonio G. Composite media and homogenization theory. New York: World Scientific; 1995; Berthier S. J Phys I 1994;4:303; Luciano R, Tamburrino A. Int J Appl Electromagn Mech 2000;11:163. [34] Bruno O. In: Bouchitte G, Buazzo G, Suquet P, editors. Calculus of variations, homogenization, and continuum mechanics. Singapore: World Scientific; 1994. [35] In the usual mean field theory due to Weiss [Weiss P. J Phys (France) 1907;6:667], a chosen atom (or spin, for magnetic problems) is viewed as interacting with the average, or mean field, of its nearest neighbors. The exact Hamiltonian is replaced by the mean field one, in which the neighbors are introduced as the mean-field acting on the central atom. However, the central atom itself influences the neighbors, and therefore the mean field usually includes a part directly attributed to the central atom. This means the mean-field includes a part that might be considered a self-interaction effect [3]. It is of some interest to note that Gubernatis and Krumhansl [Gubernatis JE, Krumhansl JA. J Appl Phys 1975;46:1875] have suggested a connection between EMT and mean-field approximations in the theory of electronic properties in disordered binary alloys. See also Mendelson KS, Schwacher D. In: Grubin HL, Hess K, Iafrate GJ, Ferry DK, editors. The physics of submicron structures. New York: American Institute of Physics; 1984. [36] Bergman DJ. Phys Rep 1978;43:377; Bergman DJ. Phys Rev B 4304;14:1976; Bergman DJ. Phys Rev B 1989;39:4598 For cases where partial information about the microstructure of the composite material is known, the spectral function can be used to derive the exact bounds on the effective permittivity. If the microstructure is exactly known, the spectral function can be used to evaluate the effective permittivity in terms of the permittivity of each constituent. A

450

[37] [38] [39] [40] [41]

[42] [43]

[44]

[45] [46] [47]

[48] [49]

[50] [51]

[52] [53]

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456 useful representation for the spectral function is obtained by expanding it in terms of power series. When the power series is outside the radius of convergence, it can be transformed into a continued fraction which converge everywhere except for the singular points.; For details on this last point, see Bergman DJ, Dunn KJ. Phys Rev B 1992;45:13262. Bergman DJ. Phys Rev Lett 1980;44:1285, Phys Rev B 1981;23:3058. Bergman DJ. Ann Phys 1982;138:78; See also Bergman DJ, Siam DJ. J Appl Math, 53, 78, and Bergman DJ. J Phys C 1982;12:4947 1979 Stroud D. Phys Rev B 3368;12:1975. Stroud D, Pan FP. Phys Rev B 1602;17:1978. Yee KS. IEEE Trans Antennas Prop 303;AP-14:1966; There have been a number of modifications to the traditional FDTD method, for solving Maxwell equations in the time domain, which avoids the limitations of Yee’s 1966 pioneering paper. One of them is implementation of the surface impedance boundary conditions to remove the body of high conductivity from the computational space and replace it with a relationship between the tangential electric and magnetic fields on the surface of the body. Since the discretization size of the space in the FDTD lattice is proportional to the wavelength of the wave propagating in the medium, removal of the medium where the wavelength becomes very small, significant simplifies the problem and decreases the requirements of the computational resource, see, for example, Maloney JG, Smith GS. IEEE Trans Antennas Prop 1992;40:38; Beggs JH, Luebbers RJ, Yee KS, Kunz KS. IEEE Trans Antennas Prop 1992;40:49. The main deficiency of the FDTD method is that only structural grids are allowed in calculation models. Consequently, it is difficult to model curved surfaces efficiently.. Taflove A, Brodwin ME. IEEE Trans Microwave Theory Tech 1970;23:623. Shlager KL, Schneider JB. IEEE Antennas Propag Mag 1995;37:39; See also Shlager KL, Schneider JB. IEEE Antennas Propag Mag 1995;37:39. The widespread use of the World Wide Web (www) offers the unprecedented opportunity of delivering information through an easy-to-use and easy-to-learn interface. Nowadays, the front-end provided by popular www browsers is known and available to millions of users, and has become a de-facto standard for the dissemination of information. Complementary and up to date information on the FDTD methods is available from the internet at http://www.fdtd.org.. Greengard L, Rokhlin V. Acta Numerica 1997;63:229; It is worth observing that using the FMM proposed in this article, to solve an electrostatic problem with the boundary element method, reduces the computational costs and the memory requirements from O(N2) to O(N), where N is the number of unknowns. See also Greengard L, Rokhlin V, Anderson C, Greengard C, Greengard L. Lecture notes in mathematics the rapid evaluation of potential fields in three dimensions. Boston: MIT Press; 1988 Liu QH. Microwave Opt Technol Lett 1997;15:158. Scaife BKP. Principles of dielectrics. Oxford: Clarendon Press; 1998. Doyle WT. J Appl Phys 1978;49:795; See also Smith ER, Tsarenko V. Mol Phys 1998;95:449 for a detailed study of the convergence of multipole expansion methods. Torquato S. Random heterogeneous materials. New York: Springer; 2001. Day AR, Grant AR, Sievers AJ, Thorpe MF. Phys Rev Lett 1978;84:2000; See also ;Roberts AP, Knacksted MA. Phys Rev E 2313;54:1996; Roberts AP, Teubner M. Phys Rev E 4141;51:1995. Brown Jr. WF. J Chem Phys 1514;23:1955; See also Brown Jr. WF. Trans Soc Rheol 1965;9:357. The first material specific (third-order upper) bound was established by Prager: Prager S. Physica, 1963;29:129. Calculation of structural parameters for third-order bounds was first discussed by Weissberg in the context of hard spheres and gradient fields from solutions to single-body electrostatic problems: Weissberg HL. J Appl Phys 1963;34:2636. Zhuck NP. Phys Rev B 1994;50:15636. Tsang L, Kong JA. Radio Sci 1981;16:303.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

451

[54] Morse PM, Feshbach H. Methods of theoretical physics. New York: McGraw-Hill; 1953; See also Romanov VG, Kabanikhin SI. Inverse problems for Maxwell’s equations. Zeist: VSP International Science Publishers; 1994. [55] Matheron G. Random sets and integral geometry. New York: Wiley; 1975; The interested reader may wish also to consult Adler RJ. The geometry of random fields. New York: Wiley; 1981. [56] Cowan G. Statistical data analysis. Oxford: Clarendon Press; 1998. [57] Tikhonov AN, Arsenin VI. Solutions of ill-posed problems. Washington: Halsted Press; 1977. [58] Kantorovich LV, Akilov GP. Functional analysis. New York: Pergamon Press; 1982. [59] A pedagogical review of this topic is presented by Baker CTH. The numerical treatment of integral equations. Oxford: Oxford University Press; 1977. [60] Savary L, Jeulin D, Thorel A. Acta Stereol 1999;18:297. [61] Jeulin D, Savary L. J Phys I France 1123;7:1997. [62] Serra J. Image analysis and mathematical morphology. London: Academic Press; 1982; See also Winkler G. Image analysis, random fields, and dynamic Monte Carlo methods. Berlin: Springer; 1995. [63] Carpenter DT, Rickman JM, Barmak K. J Appl Phys 5843;84:1998. [64] Jackson JD, Am J. Phys 2000;68:307. We would like to add that historians of science are now skeptical of the concept of discovery as a discrete event for which a particular person or persons should get credit or priority. [65] Singh ON, Lakhtakia A, editors. Electromagnetic fields in unconventional materials and structures. New York: Wiley; 2000. [66] Faraday M. Philos Trans R Soc London 1857;147:145; To get a flavor, we recommend to the interested reader Lakhtakia A. Selected papers on linear optical composite materials. Bellingham, WA: SPIE Optical Engineering Press; 1996 and the references cited therein.. [67] Maxwell JC. A treatise on electricity and magnetism, vol. 1. 2nd ed. Oxford: Clarendon Press; 1881. [68] Chew WC. Waves and fields in inhomogeneous media. New York: Van Nostrand Reinhold; 1990. [69] Stratton JA. Electromagnetic theory. New York: McGraw-Hill; 1941; See also Binns KJ, Lawreson PJ, Trowbridge CW. The analytical and numerical solution of electric and magnetic fields. New York: Wiley; 1992. [70] See, e.g., Landau L, Lifshitz EM, Pitaevskii LP. Electrodynamics of continuous media. Oxford: Butterworth-Heinemann; 1984; Jackson JD. Classical electrodynamics. 2nd ed. New York: Wiley; 1975. [71] For an entry in the numerous mixing laws which appeared in the archival literature, see Meredith RE, Tobias CW. Advances in electrochemistry and electrochemical engineering, vol. 2. New York: Interscience; 1962. [72] Van Beek LKH. Progress in dielectrics, 7. New York: Heywood Books; 1967. [73] Mossotti OF. Memorie di matematica e di fisica della Societa` Italiana delle Scienze, residente in Modena 6353;24:1850. [74] Lorenz L. Ann Phys Chem (Leipzig) 1880;11:70. [75] Lord Rayleigh (Strutt JWS) Phil Mag 1892;34:481. The analytical study undertaken by Rayleigh, in principle, would yield an exact value for the effective permittivity of the regular lattice of inclusions within the matrix, but his analytical calculation was not carried out far enough since it involves the calculation of certain static lattice sums, and his numerical study was cut short by the need of heavy computation. Following Rayleigh’s method, McPhedran and co-workers [21], and Doyle [47] computed the effective permittivity to high accuracy for the simple-body-centered-, and face-centeredcubic lattices. [76] Clausius RJE. The´orie me´canique de la chaleur, vol. 2. 2nd ed. Mons: Hector Manceaux; 1893. [77] Maxwell Garnett JC. Philos Trans R Soc London 1904;203:385, ibidem 1906;205:237. [78] Wiener O. Abh Sa¨chs Akad Wiss Leipzig Math—Naturwiss Kl 1912;32:509. [79] Lorentz HA. Theory of electrons. 2nd ed. Leipzig: Teubner; 1916.

452

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

[80] Fricke H. Phys Rev 1924;24:575; See also Fricke H. J Chem Phys 1953;57:934. [81] Lichteneker K. Physik Z 1926;27:115; See Brosseau C. J Appl Phys 1994;75:672 for a detailed discussion of the merit and an application of Lichteneker’s law of mixing in the area of impregnated polymers by nonpolar dielectric liquids; A recent paper, i.e. Zachri T, Laurent JP, Vauclin M. J Phys D: Appl Phys 1589;31:1998 has provided a theoretical basis for deriving Lichteneker’s formula; See also Goncharenko AV, Lozovski VZ, Venger EF. Opt Commun 2000;174:19. [82] Bruggeman DGA. Ann Phys Lpz 1935;24:636. [83] Biot JB, Arago F. Mem Inst France 1806;7:301. [84] Van de Hulst HC. Light scattering by small particles. New York: Wiley; 1957. [85] Milton GW. Appl Phys Lett 1980;37:300; See also Milton GW. J Appl Phys 5286;52:1981 ibid., 52, 5294, 1981; See also McPhedran RC, Milton GW. Appl Phys A 1981;26:207. [86] Ghosh K, Fuchs R. Phys Rev B 7330;44:1991. [87] Hinsen K, Felderhof BU. Phys Rev B 1992;46:12955. [88] Readers interested in the details of the phenomenology behind the GEM equation can find a comprehensive description of its development and application to composites in the following set of papers: McLachlan DS. Solid State Commun 1989;72:831; Hwang J, McLachlan DS, Mason TO. J Electroceram 1999;3:7; McLachlan DS, Hwang J, Mason TO. J Electroceram 2000;5:37 and references therein. [89] McLachlan DS, Blaszkiewics M, Newnham RE. J Am Ceram Soc 2187;73:1990. [90] Brosseau C, Queffelec P, Talbot P. J Appl Phys 4532;89:2001; See also the sequel of this paper Brosseau C. J Appl Phys 2002;90:3197 for an attempt to describe the microwave relaxation behavior of carbon black filled polymers using the generalized effective medium approach developed by McLachlan (Refs. [88,89]), which has resulted in a more general modeling capability which encompasses the limitation of Bruggeman’s analysis encountered in the vicinity of the percolation threshold. Notably good agreement between modeling and experiment was obtained in this study when the critical exponents were adjusted. [91] Sen PN, Scala C, Cohen MH. Geophysics 1981;46:781. [92] Soukoulis CM, Datta S, Economou EN. Phys Rev B 3800;49:1994. [93] Beran MJ. Nuovo Cimento 1965;38:771; See also Beran MJ. Statistical continuum theories. New York: Wiley; 1968. [94] Hashin Z, Shtrikman S. J Appl Phys 3125;33:1962. [95] Schulgasser K, Hashin Z. J Appl Phys 1976;47:424. [96] Milton GW. J Mech Phys Solids 1980;30:177. [97] Milton GW. Commun Pure Appl Math 1990;43:63. [98] Sawicz R, Golden K. J Appl Phys 7240;78:1995. [99] Corson P. J Appl Phys 3159;45:1974. [100] Broadbent SR, Hammersley JM. Proc Cam Phil Soc 1957;53:629. [101] Stauffer D, Aharony A. Introduction to percolation theory. London: Taylor and Francis; 1994. [102] Kirkpatrick S. Rev Mod Phys 1973;454:574; Kirkpatrick S. Phys Rev Lett 1971;27:1722; Clerc JP, Giraud G, Laugier JM, Luck JM. Adv Phys 1990;39:191; Isichenko MB. Rev Mod Phys 1992;64:961. [103] Sahimi M. Applications of percolation theory. London: Taylor and Francis; 1994. [104] Anderson PW. In Balian R, Maynard R, Toulouse G, editors. Ill-condensed matter. Amsterdam: North-Holland; 1979. [105] Hori M, Yonezawa F. J Phys C 1977;10:229. [106] Kawamoto H. In: Sichel EK, editor. Carbon black polymer composites. New York: Dekker; 1982 See also Miyasakara K, Watanabe K, Jojima E, Aida H, Sumita M, Ishikawa K. J Mat Sci 1610; 17:1982; Yacubowitz J, Narkis M. Polym Eng Sci 1568;26:1986.

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

453

[107] Lai JH, editor. Polymers for electronic applications. Boca Raton: CRC Press; 1994. [108] Carmona F, Barreau F, Delhaes P, Conet R. J Phys Lett 1980;41:L531; See also Carmona F, Delhaes P, Barreau F, Ordiera D, Canet R, Lafeychine L. Rev Chim Min 1981;18:498; Lagarkov AN, Sarychev AK. Phys Rev B 6319;53:1998; Lagarkov N, Matytsin SM, Rozanov KN, Sarychev AK. Physica A 1997;241:58 Observe that the rl dependence of the percolation threshold concentration in fiber filled composites was also discussed by several independent groups, raising a legitimate doubt on the general validity of the square law dependence for intermediate values of rl. Thus, there is some uncertainty over this issue.; See Boissonade J, Barreau F, Carmona F. J Phys A 2777;16:1983; Celzard A, McRae E, Deleuze C, Dufort M, Furdin G, Mareˆche´ JF. Phys Rev B 6209;53:1998; Garboczi EJ, Snyder KA, Douglas JF, Thorpe MF. Phys Rev E 1995;52:819; From the theoretical standpoint, we also also refer the interested reader to Balberg I, Anderson CH, Alexander S, Wagner N. Phys Rev B 3933;30:1984; Balberg I, Binenbaum N, Wagner N. Phys Rev Lett 1465;52:1984; Balberg I. Philos Mag B 1987;56:991; Browning L, Lodge J, Price RR, Schelleng J, Schoen PE, Zabetakis Z. J Appl Phys 6109;84:1998; Garboczi EJ, Snyder KA, Douglas JF, Thorpe MF. Phys Rev E 1995;52:819. [109] Lagarkov AN, Matytsin SM, Rozanov KN, Sarychev AK. J Appl Phys 3806;84:1998. [110] Carmona F, El Amarti A. Phys Rev B 3284;35:1987. [111] Ravier J, Houze F, Carmona F, Schneegans O, Saadaoui H. Carbon 2001;39:314. [112] Houze F, Meyer R, Schneegans O, Boyer L. Appl Phys Lett 1975;69:1996. [113] Brosseau C, Boulic F, Bourbigot C, Queffelec P, Le Mest Y, Loae¨c J, Beroual A. J Appl Phys 1997; 81:882. [114] Boulic F, Brosseau C, Le Mest Y, Loae¨c J, Carmona F. J Phys D: Appl Phys 1904;31:1998. [115] Brosseau C, Molinie´ P, Boulic F, Carmona F. J Appl Phys 8297;89:2001. [116] Tao R, Chen Z, Sheng P. Phys Rev B 2417;41:1990; See also Zhang C, Yang B, Wu X, Lu T, Zheng Y, Su W. Physica B 2000;293:16. [117] Many electromagnetic field problems can be expressed as a Freedholm integral equation of the first kind or the second kind. In integral formulations volume or surface charges or currents are usually the unknowns. These equations are solved by means of Harrington’s methods of moments, i.e. Harrington RF. Field computation by moment methods. Melbourne, Florida: Krieger Publishing Co; 1982; See also Harrington RF. IEE Proc 1967;55:136; Wang JJH. Generalized moment method in electromagnetics. New York: Wiley; 1991. It should be noted that one of the disadvantages of the Harrington’s methods of moments is the possibility of unstable solutions. Moreover, the most ill-conditioned is the Freedholm integral equation of the first kind. See also the proceedings of recent Compumags, CEFCs, and ISEMs. [118] Mur G, De Hoop AT. IEEE Trans Magn 2188;21:1985. [119] Zienkiewicz OC, Taylor RL. The finite element method. New York: McGraw-Hill; 1994. The finite element method (FEM) has considerable flexibility because arbitrary shape can be modeled. The FEM takes much more time and memory than the FDTD because it is required to solve linear equations in each time step. Consequently, the FEM is mainly used in the frequency domain and eigenmode analysis.; For a general presentation of the finite element modeling of non-destructive material evaluation, we refer the interested reader to Mackerle J. Model Simul Mater Sci Eng 1999;7:107. [120] There is a vast literature on the FEM in the context of engineering. It is also worth noting that FEM softwares with flexible mesh generators, robust equation solvers and versatile post-processors have become available in recent years. For a general introduction, the reader is refered to Strang G, Fix GJ. An analysis of the finite element method. Englewoods Cliffs: Prentice-Hall; 1973. Other important books are: Bathe KJ, Wilson EL. Numerical methods in finite element analysis. Englewoods Cliffs: Prentice-Hall; 1976, Jin J. The finite element in electromagnetics, New York: Wiley; 1993), Sylvester PP, Ferrari RL. Finite elements for electrical engineers. 2nd ed. New York: Cambridge

454

[121]

[122]

[123]

[124] [125] [126]

[127]

[128] [129] [130]

[131] [132]

[133] [134]

[135] [136]

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456 University Press; 1991, Johnson C. Numerical solution of partial differential equations by the finite element method, Cambridge: Cambridge University Press; 1987, and Hughes TJR. The finite element method. Englewood Cliffs, New Jersey: Prentice-Hall; 1987. Sto¨lzle S, Enders A, Nimtz G. J Phys I France 1991;2:401; See also Coverdale RT, Garboczi EJ, Jennings HM. Comput Mater Sci 1995;3:465; Coverdale RT, Christensen BJ, Mason TO, Jennings HM, Garboczi EJ. J Mater Sci 4984;29:1994. Weiland T. Part Accel 1984;15:245 ibidem 17, 227, 1987. This paper deals with the MAFIA (Maxwell’s equations Finite Integration Algorithm) code, which has been developed for the computation of eigenvalue equations by the Householder transformation. Bossavit A. Computational electromagnetism, variational formulations, edge elements, complementarity. New York: Academic Press; 1998; See also Bossavit A. Electromagne´tisme en vue de la mode´lisation. New York: Springer Verlag; 1993; Finlaysson BA. The method of weighted residuals and variational principles. New York: Academic; 1972 Briefly stated, the variational principle leads to an approximate solution of a partial differential equation (PDE) by extremizing a functional. The solution thus found is guaranteed to be an approximation of the unique solution of the PDE in question. Brebbia JC. The boundary element method of engineers. London: Pentech; 1980. Tuncer E. PhD thesis, Chalmers University of Technology, Gothenburg, Sweden; 2001; See also Tuncer E, Gubanski SM, Nettelblad B. J Appl Phys 8092;89:2001. Batchelor GK. Ann Rev Fluid Mech 1974;6:227; See also Schulgasser KA. Int Commun Heat Mass Transfer 1992;19:639. Insofar as the inclusions have a well-defined magnetic permeability, the problem of calculating the effective magnetic permeability is identical to that of calculating the effective permittivity. However, there is an important added aspect to the magnetic problem, i.e. the eddy current effect. Under an applied magnetic field, eddy currents are set in the inclusions which generate an induced magnetic field opposed to the applied field, and which suffer resistive losses. Thus, the eddy current effect is both diamagnetic and lossy. There are many commercial software packages available with various libraries that can be used to implement the numerical method previously described, for example, FLUX 2D and 3D (Cedrat), PHI3D, Vector Fields, Ansoft, etc. Gillespie JB, Jennings SG, Lindberg JD. Appl Opt 1978;17:989; See also Bohren CF. J Atmos Sci 1986;43:68. Calame JP, Birman A, Carmel Y, Gershon D, Levush B, Sorokin AA, Semenov VE, Dadon D, Martin LP, Rosen M. J Appl Phys 3992;80:1996. Steeman PAM, Maurer FHJ. Colloid Polym Sci 1069;270:1992; Three-phase confocal spheroid-background mixtures were also treated by Jones SB, Friedman SP. Water Resour Res 2821;36:2000. The authors would like to thank Shmulik Friedman for providing a copy of his above mentioned article. Gu G, Yu KW. J Phys D: Appl Phys 1523;30:1997. We fully agree with Kroemer [Kroemer H. Rev Mod Phys 2001;73:783] when he says that ‘‘I do not think we can realistically predict which new devices and applications may emerge, but I believe we can create an environment encouraging progress, by not asking immediately what any new science might be good for (and cutting off the funds if no answer full of fanciful promises is forthcoming’’. Frommer J. Angew Chem, Int Ed Engl 1298;31:1992. Lewis TJ. IEEE Trans Dielec Insul 1994;1:812; See also Foster KR, Schwan HP. Crit Rev Biomed Eng 1989;17:25 for a review of the dielectric properties of tissues and biological materials; Peters MJ, Stinstra JG, Hendriks M. Electromagnetics 2001;21:545 for an attempt to model the effective conductivity of human tissues, such as the cerebral cortex, the liver and blood. Polk C, Postow E. CRC handbook of biological effects of electromagnetic fields. Boca Raton, FL: CRC Press; 1986. For an overview, see Stroud D, Wood VE. J Opt Soc Am B 1989;6:778;

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

[137] [138]

[139] [140]

[141]

[142]

[143] [144] [145]

[146] [147] [148]

455

See also Stroud D, Hui PM. Phys Rev B 8719;37:1988; Zeng XC, Bergman DJ, Hui PM, Stroud D. Phys Rev B 1988;38:10970. Observe that the cubic nonlinearity is the lowest-order nonlinearity appearing in materials characterized by inversion symmetry or macroscopic isotropy. Agarwal GS, Gupta SD. Phys Rev A 5678;38:1988. Ponte Castaneda P. Philos Trans R Soc London Ser A 1992;340:531; Ponte Castaneda P, Kailsam M. Philos Trans R Soc London Ser A 793;453:1997; Ponte Castaneda P, Suquet P. Adv Appl Mech 1998;34:171; Ponte Castaneda P. J Mech Phys Solids 1991;39:45, ibidem 1996;44:827, ibidem 1997;45:317; Ponte Castaneda P, Willis JR. Philos Trans R Soc London, Ser A 1799;455:1999; The interested reader may also consult Pellegrini YP. Phys Rev B 2001;64:134211 for a self-consistent effective medium theory for random dielectric composites of arbitrary nonlinear constitutive law; Sipe JE, Boyd RW. Phys Rev A 1614;46:1992 for a theory of the enhancement of the nonlinear optical susceptibility of nano-composite materials. Yablonovitch E. Phys Rev Lett 2059;58:1987. John S. Phys Rev Lett 1987;58:2486. There exist many analogies between photonic crystals and the more familiar electronic structure of crystalline solids. In both cases, the reciprocical lattice can be described by Brillouin zones that reflect the symmetry of the real-space crystal. There is much work in this area. See, for example, Joannopoulos JD, Meade RD, Winn JN. Photonic crystals, Princeton, NJ: Princeton University Press; 1995. The periodic dielectric structure affects the dispersion relations and spatial distribution of the electromagnetic wave traveling through the photonic crystal. For large enough dielectric contrast, there exist wavelength regions (irrespective of the polarization or propagation direction) where no solutions of Maxwell’s equations in the periodic structure can be found, creating photonic band gaps. For a bibliography of papers on optical and acoustic photonic band gaps, the interested reader should refer to the URL, http://www.home.earthlink.net/jpdowling/pbgbib.html In the mid 1980s, Drexler introduced the term ‘‘nanotechnology‘‘ to describe atomically precise molecular manufacturing systems and their products. Within the past two decades, a variety of terms sharing the prefix ‘‘nano-‘‘ (from the Greek root nanos, or dwarf), such as nanoparticle, nanomaterial, nanophase, nanomachine and nanostructured, have emerged to describe certain materials, technologies, and even businesses. See Drexler KE. Proc Natl Acad Sci USA 1981;78:5275, Drexler KE. Engines of creation: the coming era of nanotechnology. New York: Anchor Press; 1986, and Drexler KE. Nanosystems: molecular machinery, manufacturing, and computation. New York: Wiley; 1992. The excitement in the field of nanotechnology has led to a sea-change in our perception of what nanodevices can do, for example, Roco MC, Williams RS, Alivisatos P, editors. Nanotechnology research directions. Dordrecht, The Netherlands: Kluwer; 2000. Most current information in this rapidly progressing area can be found on the web pages, for example, http://www.foresight.com and http://www.zyvex.com. Talbot P, Brosseau C, Konn AM. J Magn Magn Mater 2002;249:483. See also Brosseau C, Ben Youssef J, Talbot P, Konn AM. [in press]. Bichurin MI, Kornev IA, Petrov VM, Tatarenko AS, Kiliba YV, Srinivasan G. Phys Rev B 2001; 64:094409. Prinz G. Science 1660;282:1998; See also Joachim C, Gimzewski JH, Aviram A. Nature 2000;408:541. The success of ‘‘molecular electronics’’—a potential successor to microelectronics—depends on the development of so-called bottom-up manufacturing techniques, in which molecular scale components are chemically synthesized in large numbers and assembled into useful circuits. Trau M, Saville DA, Aksay IA. Science 1996;272:706. Veselago VG. Sov Phys Usp 1968;10:509. Pendry JB, Holden AJ, Stewart WJ, Youngs I. Phys Rev Lett 4773;76:1996; and the following Comment: Walser RM, Valanju AP, Valanju PM. Phys Rev Lett 2001;87:119701; See also Pendry JB, Holden AJ, Robbins DJ, Stewart WJ. J Phys Conden Matter 4785;10:1998;

456

C. Brosseau, A. Beroual / Progress in Materials Science 48 (2003) 373–456

Pendry JB, Robbins DJ, Stewart WJ. IEEE Trans Microwave Theory Tech 2075;47:1999; Smith DR, Vier DC, Padilla W, Nemat-nasser SC, Schultz S. Appl Phys Lett 1425;75:1999; Smith DR, Padilla W, Vier DC, Nemat-Nasser SC, Schultz S. PhysRev Lett 2000;84:4184. The concept of negative permeability is of particular interest because this is a regime not observed in ordinary materials, but also because such a medium can be combined with a negative permittivity to form a LHM, i.e. EH lies along the direction of k for propagating plane waves. [149] Yaghjian A. Proc IEEE 1980;68:248; See also Yaghjian A. Am J Phys 1985;53:859. [150] Bottcher CJ. Theory of electric polarization. Amsterdam: Elsevier; 1952.