Definitions of ‘effective fields’ for a randomly inhomogeneous material are offered, which guarantee automatic satisfaction of the equations of motion. The important case of a medium with periodic microstructure is included. In this special case, the definitions are completely explicit and can be applied without reference to random media. The presentation is mostly expressed in terms of electromagnetic waves. The reasoning is applicable also to other types of waves and its realization for elastodynamics is briefly outlined towards the end. Some of the effective fields are defined directly as ensemble averages, ensuring the exact satisfaction of the equations of motion, but the effective ‘kinematic’ fields to which they are related are defined more generally, as weighted averages. The main result of this work is an explicit formula for the tensor of effective properties. The important issue of uniqueness (or not) of the effective properties is explained and resolved. Self-adjointness of the original problem is not assumed. An attractive feature of the formulation is that self-adjointness at the local level implies self-adjointness at the level of the ‘effective medium’.
Metamaterials are artificially constructed composites which are designed so that they display unusual macroscopic properties, such as negative refractive index, not forbidden by physics but nevertheless not found in materials that occur naturally (e.g. Smith & Pendry 2006). The primary concern of this work is the construction of ‘effective constitutive relations’ applicable to wave propagation in such materials, but the reasoning that is employed is applicable to any composite. The major part of the presentation is for electromagnetic waves, because this reflects the majority of interest. However, a symbolic notation is employed, which is easily adaptable to other wave phenomena. Elastic waves are discussed briefly towards the end.
The definitions given for ‘effective fields’, and hence for the effective constitutive relations between them, are designed so that the field equations (Maxwell’s equations for electromagnetic waves) are satisfied without approximation. They are based on ensemble averaging, considering the microstructure of the material to display some level of randomness, but the resulting prescription can be applied to a medium with perfectly periodic microstructure, and could have been developed for such a medium without any mention of ensemble averaging. However, the formulation as given is of wider applicability.
An innovation is the introduction of weighted averages of certain fields (the electric field E in the case of electromagnetics). The use of such averages has been proposed previously (Milton & Willis 2007; Willis 2009), but here the development is more general and probably complete. It gives rise to an explicit formula for the non-local operator that defines the effective response. This is developed without assuming symmetry of the local constitutive properties, or equivalently, self-adjointness of the exact governing equations. The structure of the result demonstrates immediately that the property of self-adjointness is preserved at the level of the effective response, even when a weighted average is employed. This carries the implication, for instance, that an electromagnetic metamaterial, whose local dielectric constant and magnetic permeability tensors are symmetric, has effective response with (at most) 21 independent components rather than 36.
A significant observation, recently made by Fietz & Shvets (2010), that the effective response cannot be uniquely defined without the introduction of some ‘driving force’ additional to the usual electric current, is taken on board here. The extra driving force is needed to ensure that two fields can be generated independently. Fietz and Shvets advocated the (conceptual) introduction of magnetic monopole current; here, a magnetization tensor, whose time derivative is mathematically equivalent to (minus) the monopole current of Fietz and Shvets, is introduced instead. It has the virtue that magnetization (at least relative to a vacuum) is a well-established physical entity. In the case of elastodynamics, the quantity corresponding to magnetization is a tensor of inelastic strain, which is also a real physical entity. A brief discussion is included about the non-uniqueness that results from not admitting the additional driving force, and the further non-uniqueness that occurs if an attempt is made to infer effective response only from free waves (not driven by any current at all for electromagnetics, or driven by neither body force nor inelastic strain for elastodynamics).
In the sections that follow, first, a definition of effective response is proposed. Next, the explicit formula, equation (3.19), for the effective response is developed: this is the main result of this work. Then follows a discussion of uniqueness, after which application of the formalism to the case of periodic microstructure is explained. A brief discussion of elastodynamics is given next. The paper is concluded with a few general comments. The symbolic notation that is employed has the virtue of displaying the logical structure, but at the expense of covering some of the detail. This is compensated by including in appendix A derivations of some of the key formulae, employing explicit suffix notation.
2. Definition of effective response
Maxwell’s field equations are adopted in the form 2.1 The problem that will be considered is for a medium occupying a domain Ω, with linear constitutive relations 2.2 Here, ε and μ are second-order tensor convolution operators with respect to time, which depend on position x. The asterisk denotes the operation of convolution with respect to t. In this sense, ε and μ are functions of x and t. In the frequency domain, they become functions of x and radian frequency ω. The addition of the magnetization M is unconventional for two reasons. First, M is taken relative to the medium under consideration, not to a vacuum. Second, it will be assumed that M is imposed by some external means, in the same way that the current j in equation (2.1)2 is the externally imposed current. The reason for its inclusion will be explained later.1
The problem is specified by the requirements that equations (2.1) and (2.2) must apply in Ω. In addition, initial conditions are specified throughout Ω (the simplest are that all fields are zero for all t<0) and a suitable set of conditions are applied on the boundary ∂Ω of Ω. Although more general conditions could be specified, it is sufficient for the present purpose to divide ∂Ω into complementary parts SE and SH. Over SE, tangential components of E are given, while tangential components of H are specified over SH. Thus, 2.3 while 2.4 where n denotes the outward normal to the boundary.
The first Maxwell equation (2.1)1 is satisfied identically if a vector potential A is introduced, in terms of which 2.5 The problem then reduces to that of finding A(x,t).
Now to define effective relations, the medium is taken to be random. Thus, the constitutive functions become functions of x, t and a parameter α that belongs to a sample space over which a probability measure p is defined. It is perhaps relevant to remark straightaway that a medium with periodic structure, with period cell Q say, fits within this class if it is recognized that the exact location of a specified point in one cell is not known exactly but is placed at random. In this special case, the parameter α can be identified with the spatial coordinate of that point, and it is uniformly distributed over Q. It is desirable that Maxwell’s field equations should be satisfied exactly by the effective fields, Eeff, etc. Satisfaction of the first equation, (2.1)1, is ensured by requiring that 2.6 for appropriate choice of Aeff (discussed below). Taking the ensemble mean of the second Maxwell equation, (2.1)2, yields 2.7 where the ensemble mean of a random function f is defined as 2.8 Effective relations are then defined by expressing Heff=〈H〉 and Deff=〈D〉 in terms of Aeff, and hence in terms of Eeff and Beff.
The simplest and most obvious prescription for completing this programme is to define 2.9 and associated with this, to take the current j, the magnetization M and the functions E0 and H0 that define the boundary data to be sure, i.e. the same for every realization of the medium. It is possible, however, to develop a more general prescription, which includes equation (2.9) as a special case, by adopting for Aeff a weighted mean 2.10 where w(x,α) is a random field with 〈w〉=1. Considering the relation (2.6)2, the choice (2.10) is equivalent to adopting the weighted mean 〈wE〉 of E. This generalization is more than just academic: for a metamaterial comprising a uniform matrix in which is embedded an array of complex resonators, it may well be appropriate physically to define Eeff as an average over just the matrix, in which case w would be uniform over the matrix but zero over the regions occupied by the resonators. This prescription was developed via rigorous analysis of the ‘homogenization limit’ (or long-wavelength limit) for a particular problem of scattering by a slab composed of material with periodic microstructure by Kohn & Shipman (2008). Here, following Willis (2009) (see also Milton & Willis 2007), it is proposed for a general random medium, and for any wavelength or frequency. The section that follows develops effective constitutive relations in terms of the weighted mean vector potential , including as a special case the unweighted mean, which is obtained by setting w(x,α)=1.
3. Formal derivation of effective relations
It is convenient to consider the problem in Laplace-transformed form, so that ∂/∂t is replaced by s and convolutions with respect to time become ordinary products. Dependence on the transform variable s will not be shown explicitly. The structure of the reasoning to follow is applicable not only to electromagnetic waves but also to other waves, including elastodynamic waves. The structure is highlighted by introducing a symbolic notation as follows: 3.1
It follows that the partial differential equation governing A can be written as 3.4
The vector potential A so defined can be expressed, at least formally, in terms of Green’s function G associated with equation (3.4) and the boundary conditions. This is a 3×3 matrix with components Gij(x,x′) whose adjoint G† (with ij component ) satisfies the equation2 3.6 Here, I is the 3×3 identity matrix, δ(x−x′) is the three-dimensional Dirac delta, and the curl operation within operates on the first argument that it encounters. Homogeneous boundary conditions corresponding to equation (3.5) are applied to G†. The required representation is 3.7 An integration by parts reduces this to the form 3.8 having taken account of the boundary conditions for both A and G†. The exact meaning of the representation (3.8) is clarified in appendix A, where it is derived explicitly, employing suffix notation.
In correspondence with employing the constraint that 〈wA〉 is specified, it is desirable to specialize the current j and the field H0 that defines the tangential component of H over SH to have the form 3.9 where j1 and H1 are sure. Thus, equivalently, 3.10 However, the magnetization M (equivalently, m) and the function n×A0 that define the boundary condition over SE are taken to be sure.
In the sequel, equation (3.8) will be written in the symbolic form 3.11 where the definition of A0(x) is extended from SE to the whole of Ω by the formula 3.12 It is required now to calculate 〈wA〉. From equation (3.11), 3.13 having exploited the form (3.10) of j and H0. It is possible to express A0 as given by equation (3.12) in an important alternative form. Since n×A0 is sure on SE, it can be replaced by n×〈wA〉 there. Furthermore, having made this replacement, the integral over SE can be extended to an integral over the whole of ∂Ω, on account of the homogeneous boundary condition satisfied by G† on SH. Transforming this boundary integral to an integral over Ω then gives 3.14 since G† satisfies equation (3.6). Equation (3.14) is also derived, employing suffix notation, in appendix A. It follows, since m is sure, that 3.15 and 3.16 Hence, reverting to equation (3.11), 3.17 Now, and, therefore, 3.18 where the effective constitutive operator is given by 3.19 It is evident that equation (3.19) defines as spatially non-local; inverting the Laplace transform demonstrates that is also a convolution operator with respect to time, even in the special case that the local operator is independent of s (or t in the time domain).
The effective constitutive relation (3.18) has the general bi-anisotropic form 3.20 The constitutive matrix is not self-adjoint in the general case of non-symmetric , but when is symmetric (i.e. ε and μ are symmetric), the problem is self-adjoint, G†=G and , 3.21
4. Remarks on uniqueness
It was admitted from the outset that the magnetization M was introduced artificially. Its purpose is clear from the structure of the effective relation (3.18): it ensures independence of the two terms in the ‘column vector’ so that no alternative relation is possible. The proposal of Fietz & Shvets (2010) to introduce into the first Maxwell equation a magnetic monopole current has the same effect. The reasons for the present advocacy of introducing magnetization M are that (i) it is at least a familiar physical entity relative to a vacuum and corresponds to a distribution of microscopic loops of electric current and (ii) there is an analogous quantity, called ‘eigenstrain’ or ‘inelastic strain’ in elasticity (discussed in §6).
The consequences of not admitting M (or magnetic monopole current) are now briefly discussed. Suppose not only that M≡0, but also that the boundary condition on SE is homogeneous so that the integrals involving n×A0 over SE disappear.3 The development of the preceding section can still be followed unchanged, to give the expression (3.19) for . However, it is possible to add any operator, say, for which 4.1 for all 〈wA〉 compatible with its boundary condition. Spelling this out, 4.2 Thus, 4.3
It is perhaps relevant to remark that Willis (2009) developed the notion of effective relations using weighted mean fields, but did not introduce any field corresponding to M, nor n×A0 over SE. In the present notation, this has the effect of reducing equation (3.13) to 4.4 while equation (3.11) becomes A=Gw〈j〉. It follows that 4.5 For a one-dimensional example (which is automatically self-adjoint), Willis (2009) manipulated this into a relation of the form (3.18) (with m=0), which was also self-adjoint. However, as shown above, this form is not determined uniquely by equation (4.5).
Suppose now that not even electric current j is admitted. It is then only necessary that equations (4.1) should apply, for all 〈wA〉 that satisfy the equation of motion 4.6 Thus, it is possible to deduce solely from observations of free waves4 that, if is an operator for which for all 〈wA〉 that satisfy equation (4.6), then so does , where 4.7 for an arbitrary operator Λ, where has the form given in equation (4.3). It may be noted that has the form of a row vector, and Λ has the form of a column vector. In this sense, the term has the form of a rank-one matrix.
5. Implementation for an infinite periodic medium
It has already been remarked that a periodic medium can be regarded as random if the exact location of any one period cell is considered to be unknown. To express this precisely, suppose that a function F0(x) is periodic with period cell Q. Then, make the definition 5.1 The class of random fields of concern here are functions Fy(x), where y plays the role of the variable α and is uniformly distributed over Q.
Thus, the random medium is taken to have constitutive tensor 5.2 where is Q-periodic and y is uniformly distributed over Q. The weight w(x) will be similarly defined. The task at hand is to evaluate the expression (3.19) in this special case. This requires the calculation of Green’s function G. For the realization y=0, call Green’s function G0(x,x′). Then, for general realization y, 5.3 The required ensemble averaging reduces to integration with respect to y over the period cell Q. Explicitly, equation (3.19) gives 5.4 where 5.5 5.6 and 5.7 How next to proceed depends on the form in which Green’s function G0(x,x′) has been obtained. One way is to express it as a superposition of Floquet–Bloch waves, but this is not essential.
Willis (2009) proposed, and implemented this scheme for a one-dimensional example, but in the absence of a field M(x). The ensemble averaging was done on equation (4.5). The effective tensor so derived was not unique, though this was not noted at the time.
There is, however, an alternative procedure for a periodic medium, which is equivalent. The effective tensor of a medium of the type under consideration is translation invariant, and can be investigated via its Fourier transform 5.8 Equivalently, it suffices to consider the problem, defined in realization y by equation (3.4), with 5.9 where and are constants. In this case, equation (3.4) supports a solution of the form5 5.10 where the solution for realization y=0 is A0(x)e−ik⋅x. Hence, h(x) has the similar form 5.11 It follows immediately that 5.12 where 5.13 is the mean value of the periodic function h0. It is also obvious, independent of the notion of ensemble averaging, that 5.14 where , because for any y. The calculation of 〈wA〉 is similar, 5.15 These definitions give and , and hence also , as functions of and , whose elimination yields the (unique) effective relation , where .
It may be noted that this prescription for defining effective properties only requires the solution of the specified problem in one single realization of the medium. It was presented in the restricted context of free waves by Nemat-Nasser et al. (in press). Fietz & Shvets (2010) proposed the solution of equations driven by a current (in the present notation) and a similar magnetic monopole current. Their innovation was that this enabled the deduction of a unique effective property operator, . They emphasized that their approach could be followed in conjunction with any averaging scheme, but the one that they employed was based on volume averaging of the full fields rather than just their periodic parts, and included dividing some of the integrals by an ‘effective volume’, to avoid obvious violation of the Maxwell equations for the effective fields; see also Smith & Pendry (2006). By contrast, the present averaging prescription leads to exact satisfaction of Maxwell’s equations. Explicit examples demonstrating this in practice, in the context of free waves, are given by Nemat-Nasser et al. (in press).6
In the case of elastodynamics, the equation of motion (which is analogous to the second Maxwell equation (2.1)2), Laplace transformed with respect to time, is 6.1 where σ denotes the Cauchy stress tensor (with components σij), p denotes momentum density (a vector, with components pi) and f (with components fi) is body force. In terms of the displacement vector u (which is analogous to the vector potential A), the constitutive relations can be written 6.2 where C denotes the fourth-order tensor of moduli and ρ is the mass density. The eigenstrain η has been introduced here, for the same reason that the magnetization M was introduced in the electromagnetic case. It can be realized in practice (with limitations) by thermal distortion, or plastic deformation, or phase transformation. The tensor C displays the symmetries 6.3 so that stress depends on displacement gradient ∇u only through the strain 6.4 It also suffices to take η to be symmetric. The tensor C could depend on s so that the body is viscoelastic. It may also display the symmetry 6.5 but this is not assumed. The problem will be for a body occupying domain Ω, with displacement u=u0 specified on a part Su of the boundary ∂Ω, and traction n⋅σ=n⋅σ0 on the complementary part ST. The initial conditions for the corresponding problem in the time domain are that all fields are zero for t<0.
The entire problem can be placed in the form of equation (3.4) by making the following definitions: 6.6 where I is the 3×3 identity. All of the reasoning that has been developed applies equally when these relations are in force.7 In particular, the effective constitutive operator is given by equation (3.19). The form of the constitutive relations to which it leads may be expressed as 6.7 as proposed by Willis (2009). In contrast to the original density ρ(x), the ‘effective density’ is a second-order tensor-valued non-local operator. This was demonstrated explicitly in the case w=1 by Willis (1985). Effective relations of the form (6.7) (with w=1) were developed by Willis (1997), but they were implicit in much earlier work, e.g. Willis (1981). A first attempt to allow for a weighted average of displacement was made by Milton & Willis (2007). When the tensor of moduli C has the symmetry (6.5) necessary for self-adjointness, the effective relations (6.7) are also self-adjoint.
7. Concluding remarks
A brief comment on the relevance of any scheme based on ensemble averaging may be in order. The ensemble average, whether weighted or not, of any field does not represent directly any local average of that field. Suppose, however, that the equations governing the weighted mean field 〈wA〉 have been solved. Equation (3.17) expresses the exact field A in any realization, in terms of 〈wA〉. The formula may be hard to use because it involves Green’s function of that one realization, but at least equation (3.17) shows that the local field can, in principle, be reconstructed, if the mean field 〈wA〉 is known. If the medium is periodic, the situation is better because Green’s function is known for any realization if it is known for one; see equation (5.3). For more general random media, some approximation scheme will be required. Several such schemes do exist, usually based on a formulation relative to a uniform comparison medium.8 See, for instance, Tsang & Kong (2001). If G0 is Green’s function for a body with constitutive tensor , the relevant exact formula is 7.1
The present theory is exact, in the sense that not only Maxwell’s equations but also the boundary conditions are satisfied, because Green’s function G is that for the original boundary-value problem. It remains likely in practice, however, that calculations will be performed based on the corresponding Green function for an infinite medium. The resulting constitutive relations are then not applicable in some boundary layer whose thickness is of the order of the microscopic length scale. If attention is restricted to problems in which the generated mean wave is composed predominantly from wavelengths significantly greater than the microscopic scale, the ‘infinite-body’ constitutive operator closely approximates a local operator that can formally be applied right up to the boundary. The error so induced is small and tends to zero as the ratio of microscale to predominant wavelength tends to zero. This is the homogenization limit analysed by Kohn & Shipman (2008). For smaller wavelengths, however, this aspect merits further attention, not just in the context of the dynamics of metamaterials; it is beyond the scope of the present work.
The main result of this paper is equation (3.19). It is applicable not only to electrodynamics, but also to other wave-propagation phenomena including elastodynamics. It results from a definition of effective response that leads to exact satisfaction of the equation of motion (Maxwell’s second equation for electrodynamics). This is because the effective fields that appear in this equation are ensemble averages of the corresponding local fields in each separate realization. These effective fields are related constitutively to fields obtained from a weighted average of the corresponding local field (vector potential A for electrodynamics or displacement vector u for elastodynamics), thereby satisfying automatically the ‘compatibility relations’ (Maxwell’s first equation for electrodynamics). The prescription of course includes the possibility of employing an unweighted average for A or u, by taking the weight w(x)≡1. Weighted averages may well be useful in practice for metamaterials, for which it may be desirable to perform averaging, excluding regions containing resonating inclusions. Ensemble averaging may perhaps be viewed as somewhat esoteric, considering that most metamaterial designs are likely to feature periodic microstructure. However, it has been shown how periodic media can be subjected to the proposed ensemble averaging by exploiting the Floquet–Bloch structure of the wave fields, meaning that the fields have only to be constructed (or computed) in a single realization of the medium. It has further been shown how the effective relations can be obtained by averaging just the periodic parts of the relevant fields in this case, with no explicit reference to ensemble averaging.
The formulation picked up on a recent observation by Fietz & Shvets (2010) that effective relations can be uniquely defined only if some additional driving force is introduced. In the electromagnetic context, Fietz and Shvets proposed the introduction of magnetic monopole current. The extra field introduced here is a field of magnetization for electromagnetics (or a field of eigenstrain or inelastic strain, for elasticity); it serves the same purpose and preserves Maxwell’s first equation in its usual form (2.1)1.
The proposal of Fietz & Shvets (2010) was to calculate the averages of the participating fields, and then solve a suitable set of linear equations to deduce the coefficients that relate them. The explicit formula (3.19) is believed to be new, even in the case w=1.
One immediate implication of equation (3.19) is that, if the original problem is self-adjoint, then the effective medium formulation is also self-adjoint.9 In the electromagnetic context, this results in there being 21 effective constants rather than 36. A further implication is that, if the original problem is self-adjoint and so permits a variational formulation, there is a corresponding variational formulation for the effective medium. Full discussion of these variational formulations will be presented elsewhere.
Appendix A. Some derivations in explicit suffix notation
The governing equation (3.4) for A is A1 where eijk is the alternating symbol and ∂j means ∂/∂xj. Equation (3.6) is A2 Multiplying equation (A1) by and equation (A2) by An, and subtracting one from the other gives A3 Integration over Ω and use of the divergence theorem then gives A4 since Hk=(μ−1)kl(elmn∂mAn−Ml). Equation (A4) is equivalent to equation (3.8), taking account of the boundary conditions satisfied by G† and the fact that .
In the context of elastodynamics, the realization of equation (3.4) is A6 and equation (3.6) governing G† becomes A7 Multiplying equation (A6) by and equation (A7) by uk and subtracting one from the other gives A8 Integrating over Ω and employing the divergence theorem gives A9 This provides the analogue of equation (3.8) for elastodynamics, taking into account the boundary conditions for G†. The derivation of equation (3.14) in the context of elastodynamics mirrors the reasoning that lead to equation (A5) above, and is not repeated.
↵1 Relative to the present notation, Fietz & Shvets (2010) adopted the constitutive relation B=μ*H, in conjunction with the Maxwell equation , where j(m) is the magnetic monopole current. Thus, j(m) is mathematically equivalent to .
↵2 With G defined via its adjoint in this way, the reasoning that follows requires no assumption of symmetry of the constitutive tensors ε and μ.
↵3 This is closely related to disallowing a magnetic monopole current on the surface, in the same way that specifying n×H is related to surface electric current.
↵4 In practice, ‘observations’ are likely to be the results of detailed computations.
↵6 The exact context was one-dimensional elastodynamics, but this is isomorphic to the corresponding problem in electrodynamics.
↵8 In fact, preliminary work employed a comparison medium, until it was realized that it was equivalent to direct use of the actual medium and its Green’s function G, for which the representation for A took the relatively concise form (3.17), leading to the final result (3.19).
↵9 Of course, this depends on adoption of the present definition of effective response.
- Received November 29, 2010.
- Accepted January 10, 2011.
- This journal is © 2011 The Royal Society