Homogenization of periodic metamaterials by field averaging over unit cell boundaries: use and limitations

We revise the method of periodic metamaterials homogenization initially proposed by Pendry, Holden, Robbins and Stewart (PHRS). The shortcomings of the PHRS derivation of the basic formulae of their method are outlined, subtleties of the method implementation are discussed and the range of validity of both the PHRS method and its later modifications are analyzed. We then give a rigorous derivation of the PHRS averaging formulae in the static approximation and modify the PHRS method to account for the phase advance of an incident wave across the unit cells of metamaterials beyond the quasistatic regime. The advantages of our proposed method are illustrated by numerical calculations of the effective parameters of some periodic metamaterials.


Introduction
Metamaterials for electromagnetic and optical applications are artificial composites typically consisting of subwavelength inclusions (e.g. spheres, rods, swiss rolls or split ring resonators (SRRs)) arranged periodically in a host dielectric medium. On a scale much larger than the period of the metamaterial's lattice, they can be treated as effective homogeneous media characterized by their effective material parameters-effective permittivity ε eff , permeability µ eff and conductivity σ eff .
Calculation of the effective parameters of a metamaterial from its structure is of great importance in metamaterials design and has been accomplished using various homogenization theories and approaches that have been proposed previously (see e.g. [1][2][3][4][5][6][7][8][9][10][11][12]). For example, the effective permittivity of a diluted lattice of particles in the quasistatic regime (when the size a of the lattice's unit cell is small compared to the operating wavelength, a λ) can be calculated by using a simple Maxwell Garnett formula [2]. However, beyond the quasistatic regime and for high volume fill factors more sophisticated approaches must be employed.
Advances in computational electrodynamics have enabled numerical simulation of electromagnetic wave propagation through complex structures such as metamaterials and in turn have allowed the effective parameters of these structures to be calculated directly from simulation. One popular method of calculating effective material parameters from simulated data is the transmission line equivalent model (TLEM) [13,14]. However, implementation of this approach is at best very difficult in the case of metamaterials that are only one to five layers thick in the direction of wave propagation and whose lattice constant is in the range of 0.1λ 0 -0.4λ 0 (see [15] and references therein). However, it is on these materials that contemporary work on metamaterials has often focused.
Another particularly appealing method of numerical homogenization of metamaterials has been proposed by Pendry, Holden, Robbins and Stewart (PHRS) [6] (a detailed review of the method is given in [16]). One of the key advantages of this method is that it calculates the effective permittivity and permeability of periodic metamaterials with simple cubic (SC) lattices by using only information on local E, D, H and B fields on the boundaries of a unit cell of a metamaterial. The local field values can be obtained by numerically solving the corresponding boundary-value problem for Maxwell's equations with an appropriate solver.
When implementing the PHRS method, the average values of the fields appearing in the definitions of the effective permittivity and permeability are defined as either surface integrals of the local fields over the faces of a unit cell (for D and B fields) or as line integrals of the local fields over the edges of the cell (for Eand H fields). The average value of each Cartesian component of each field is calculated by using only one face (for D and B fields) or one edge (for E and H ).
Although the PHRS method has been claimed to be valid beyond the quasistatic regime and for virtually any contents of the unit cell [16], it returns, in fact, implausible results even in simplest cases of lossless lattices: it returns complex valued and increasingly inaccurate material parameters as the size of the unit cell increases [17,18]. Also, this method suffers from not being properly justified theoretically: the adopted definitions of the field averages are rather arbitrary and the statements regarding the range of validity of the method are contradictory (a more detailed discussion is given in section 2.1).
Later modifications of the PHRS method allow one to obtain, in some cases, more plausible results either by removing a phase factor from the calculation results [16,17] or by averaging the local fields simultaneously over several faces (or edges) of a unit cell [9]. Still, these modifications completely depend on the PHRS averaging procedure and do not provide a theoretical validation of the procedure itself. Also, the results from the different modifications do not always comply with each other and with the corresponding results from other homogenization theories 4 .
Yet another problem not analyzed previously is the relationship between the conventional (volumetric) field averaging and the averages used in the PHRS method. The former was first introduced by Lorentz [19,20] and then modified, by invoking smoothing functions, or mollifiers, by Russakoff [21] (see also [22][23][24]). Spatial averaging has been widely adopted and has provided fruitful results in various branches of physics such as macroscopic electrodynamics, optics, solid state physics and effective media theories [19,[22][23][24][25][26][27][28][29][30]. Regardless of these successes, volumetric averaging has been considered by some researchers as inapplicable in the case of metamaterials made of thin wires or sheets of metal [6] or even for metamaterials at all [9,31]. Interestingly, this apparent inapplicability motivated PHRS to introduce their own definitions of the field averages. Since a detailed explanation was not provided, it is still not clear why volumetric averaging is not valid in the case of metamaterials. Lastly, the authors of the PHRS method did not consider whether (and in which way) their method could be obtained from conventional volumetric averaging method.
The aim of this paper is to review the PHRS averaging method, provide a rigorous derivation of its basic formulae from first principles, discuss subtleties of its implementation, modify it to extend its range of validity and analyze its limitations.
Throughout the paper, we consider an infinite medium consisting of uncharged inclusions of arbitrary shape and material located in air or vacuum. The inclusions are arranged in an SC lattice with the unit cell size a. As in [6], the inclusions do not intersect the boundaries of the unit cell. Also, it is assumed that there are no impressed currents in the lattice (J ext = 0). This paper is organized as follows. In section 2, we analyze the drawbacks of the PHRS derivation of the field averaging formulae for lattices of particles and present a rigorous derivation of the formulae in the static approximation. Then, we try to modify the PHRS formulae to make them more applicable to lattices with appreciable unit cell size. In section 3, numerical results and a discussion are given to illustrate the accuracy and the range of validity of both the initial and modified methods of field averaging.

Theory
In this section, we first outline the averaging method proposed by PHRS and then give, within the static approximation, its theoretical justification based on conventional volumetric averaging.
In general, the tensors of the effective relative permittivity ε eff and permeability µ eff of inhomogeneous media are defined via constitutive relations [32] Here, ε eff ≡ ε 0 ε eff,r , µ eff ≡ µ 0 µ eff,r (ε 0 , µ 0 are permittivity and permeability of free space and ε eff,r , µ eff,r are relative effective permittivity and permeability of medium) and the angle brackets · · · denote the spatial (volumetric) averages of the respective quantities, V being the volume of averaging. For periodic structures-e.g. lattices of particles-V is merely the volume of the lattice's unit cell, V = V cell [28,29]. Note that the quantities E , D , H and B appearing in equations (1) are the average values of the respective macroscopic fields E, D, H and B (see [9,16]). The averaging (2) is performed on the scale of the lattice constant (say ∼1 cm for metamaterials operating in GHz region) and is different from what is traditionally used to obtain macroscopic Maxwell's equations from their microscopic counterparts-Maxwell-Lorentz equations-and to introduce the macroscopic fields themselves. In the latter case, averaging over a 'physically infinitesimal' volume (e.g. a unit cell of a crystal), typically a few nanometers in size, is assumed. Also, the averaging (2) is the simplest possible spatial averaging: it does not involve any smoothing functions (see [21]).
From equations (1) it follows the coordinate representation of ε eff and µ eff , with F i being the average ith component of the respective field F (F = E, D, H, B). When considering the effective parameters of periodic structures made of thin wires or sheets of metal, PHRS noted [6] that calculations based on (1)-(2) should always give ε eff,r = µ eff,r = 1. (This issue will be resolved in section 2.2.5.) They then questioned the definitions of the average values of E, D, H and B fields one has to use in calculation of ε eff and µ eff according to (3). The integration path L z and surface S z to calculate the average values of z-components of the fields according to (6)- (9). To calculate the averages of ith component (i = x, y, z), the line L i and the normal n i of the surface S i are to be along ith-axis.

Pendry, Holden, Robbins and Stewart (PHRS) method of field averaging
The PHRS averaging formulae are obtained [6,16] from the integral form of Maxwell's curl equations applied to the faces and the respective bounding edges of a unit cell, see figure 1(a). By interpreting the line and surface integrals as quantities that define the average values of the respective fields, PHRS introduced the averaging formulae where the integration path L i and surface S i (i = x, y, z) are chosen as shown in figure 1(b). It might seem that the PHRS definitions (6)- (9) are an obvious consequence of the exact Maxwell's equations. A closer examination, however, reveals that this is not true. In fact, the definitions are a result of an arbitrary interpretation of the line and surface integrals appearing in (4) and (5). Also, the approach used by PHRS to derive the averaging formulae suffers from drawbacks which are analyzed below. Firstly, the above definitions are drawn from Maxwell's curl equations (4) and (5) applied to a single face of a unit cell. Yet, it is still not convincing why such defined field averages have to satisfy not only the equations they are drawn from but also the constitutive relations (1) and, thus, define the effective constitutive parameters (ε eff , µ eff ) of the lattice.
Secondly, definitions (6) and (7) are drawn from separate terms of the contour integrals appearing in the left-hand sides of equations (4) and (5), respectively, while in definitions (8) and (9) the entire surface integrals from the right-hand sides of the same equations are used. This leads to a rather unexpected result: the field averages do not satisfy Maxwell's equations from which they were derived. For instance, in case of S chosen as in figure 1(a), equation (4) reads The four line integrals in the left-hand side define four PHRS averages of E field over respective edges L 1 , . . . , L 4 . It is only the sum of these four that equates to the right-hand side; none of the line averages individually satisfies (4). This inequality contrasts sharply with that of Lorentz volumetric averaging where the averaged fields do satisfy macroscopic Maxwell's equations.
Thirdly, the choice of the integration path L or surface S is not unique: generally speaking, one might choose between four edges or two faces of the unit cell to calculate the average of a specific field component. However, PHRS did not provide any reason for their choice illustrated in figure 1(b) nor did they argue that all possible choices are equivalent. Moreover, in a later work [16] Smith and Pendry employed another choice of L which does not comply with the PHRS prescription (this will be discussed in detail in section 3). Recently, it has also been proposed [9,31] to use all 12 edges of the unit cell in calculation of the averages of the 'curlconforming' fields E, H and all the six faces of the cell for the 'div-conforming' fields D, B. Hence, it is not clear whether the averages over different edges or faces are equivalent to each other and, if not, which of the averages do represent the true average values of the fields and are to be used in calculation of ε eff , µ eff .
Fourthly, one may question whether or not the PHRS definitions (6)- (9) are related in any way to the conventional definition (2) of the field averages. The way in which the definitions (6)- (9) have been introduced by PHRS does not provide insight as to this possible relation. Further, PHRS noted [6] that for lattices of thin wires or sheets of metal, the conventional definition would always give relative values ε eff,r = µ eff,r = 1, the same as for free space. In contrast, their method of field averaging is reported [6,16,17] to give plausible results for the lattices of this type. Therefore, one might conclude that the PHRS averaging is not related in any way to the volumetric one and is even preferable to the latter.
Finally, the PHRS approach does not allow one to define the range of validity of their averaging formulae: On the one hand, the integral Maxwell's equations which were used to introduce the field averages are valid for any size of the integration areas and contours. The finite-difference equations replacing, after applying the PHRS averaging procedure, the integral Maxwell's equations (see (8) in [16]) are exact in the same sense and to the same extent as Maxwell's equations themselves, since the former are nothing but an abbreviated form of the latter. Respectively, the PHRS approach and its modifications were claimed to be correct for any unit cell size [16,31]. On the other hand, such claims contradict to the general result that the very material homogenization concept is valid for only the case of a λ [28,29,33].
Overall, the approach used by PHRS to derive the averaging formulae seems to be misleading although the results they obtained may still be correct, at least in some cases. For example, averaging Maxwell's curl equation ∇ × H = ∂D/∂t by applying, according to PHRS, the line averaging · · · L to its left-hand side and the surface averaging · · · S to the right-hand side, results in an equation for the average quantities H L and D S that is not true, in general case. However, in the static case the averaged equation is still true since ∇ × H = ∂D/∂t = 0 for static fields in the absence of impressed currents. Based on this observation, below we try to rigorously validate the PHRS averaging formulae (6)- (9) in the static approximation starting from the well-established conventional definition (2) of the volumetric average (which works fairly well in the long-wavelength regime) and by using the general properties of static electric and magnetic fields.

Derivation of PHRS averaging formulae
Within our approach, the problem of validation of (6)-(9) can be formulated schematically as and is thus reduced to validating the following two sets of replacements: 1. replacement of the volumetric average by a surface one, · · · V → · · · S , for D and B fields and 2. replacement of the volumetric average by a line one, · · · V → · · · L , for E and H fields.
Below we shall make use of the following basic properties of static fields in a periodic structure: 1. The properties of static electric and magnetic fields as described by Maxwell's integral equations where C and Sare arbitrary but closed path and surface respectively, Q is the total charge enclosed by S. 2. Periodicity of total fields in an infinite periodic structure where F denotes any of the four fields and a = n 1 a 1 + n 2 a 2 + n 3 a 3 is the translation vector (a i primitive translation vectors of the lattice, n i integers, i = x, y, z). Since we consider only SC lattices, the norms of all the vectors a i will be set to a i = a.
Let us first consider the replacements for E and H fields. (10) the latter is due to the field periodicity (14). After changing the direction of integration in the second integral in (15) and, simultaneously, the sign before the integral, equation (15) reads Equation (16) tells us that the line integral of the E field along a line segment parallel to the z-axis does not depend on the position of the segment within the unit cell. In other words, the line integral between two corresponding points A and D has the same value for any choice of the points: We may, therefore, use any such segment-lying either on an edge, a face or inside the cell, see figure 2 The coordinate form of (17) can be obtained by noting that E · dl = E l dl (only the longitudinal component E l of E along L z contributes to the integral). Thus, we have Since there is nothing special about the z-direction in (10), by considering closed contours with their line segments oriented in either the xor y-directions, we can obtain where L i denotes line segments parallel to ith-axis (i = x, y, z) and C i are constants. Now we can use the constancy of the line integrals to reduce the volumetric averages E i V of E i components to their line averages E i L . Reducing the triple integral appearing in the and accounting for (19) and V = a 3 , we may rewrite E i V as According to (20), volumetric average of the ith component of a static and periodic E field over the volume of the unit cell equals to the line average of that component calculated over a line segment L i of length a parallel to the ith-axis. The arbitrariness of the location of L i means that PHRS requirement that L i must be one of the edges of a unit cell is excessive.
The current density J in (11) consists, generally speaking, of two terms J ext and J ind are the densities of external (impressed) and induced currents. The former is assumed to be zero, J ext = 0, see section 1. The latter is related to the local value of the total electric field through Ohm's law σ is the local value of the conductivity of the inclusion material. In the static case considered here, electric field E cannot induce steady currents flowing through inclusions of finite size. Thus, J ind = 0, J = 0, and (11) now reads which is similar to (10) for E field. Then, a consideration completely similar to that given in the previous section for the E field (or formal replacement E → H in (20)) gives (13) (13) as In (23), dS ≡ n dS with n being the unit vector of the external (relative to the enclosed volume) normal to the surface element dS. For example, on the surfaces S bott and S mid dS bott ≡ n bott dS z = −ẑ dx dy, dS mid ≡ n mid dS z =ẑ dx dy withẑ being the unit vector in the z-direction, and dS z ≡ dx dy is a surface element normal to the z-axis.
Due to the periodicity of the B field in a lattice, the field vectors at corresponding points on opposite faces of a unit cell are equal (see figure 3), which leads to or, in a coordinate form Accounting for the arbitrariness of S mid , this gives with S z being the arbitrary cross section of the unit cell whose normal is in the z-direction.
Following an identical procedure (or formally replacing the z index in (25) by x or y), one can obtain two more equations that can be written, along with (25), in a general form the constants being, generally speaking, different for different i values (i = x, y, z). The constancy (26) of the surface integrals of the B field can now be used to reduce the volumetric averages B i V of the B components to their surface averages B i S , the procedure is analogous to that used in the derivation of (20). In this way, we obtain with S i not necessarily located on a face of the unit cell, in contrast to what is required within the PHRS approach.

Replacement D i V → D i S .
Care needs to be exercised while replacing D i V by D i S . Proceeding in full analogy with the above derivation for B field, we could start from Gauss's theorem for D field written as (29) or, accounting for zero contribution of the side fluxes D side, j due to D field periodicity,  (34) and (36). Volume V (z) contains total induced charge Q(z) in it; a layer of thickness dz next to the surface S z (z) possesses total charge dQ(z).
where Q free is the total free charge enclosed by the surface S. Unlike the respective equation (24) for the B field, (30) has a non-zero right-hand side. Therefore, we cannot conclude that D flux in the z-direction is always conserved within the unit cell. For lattices of dielectric inclusions, D flux is indeed conserved due to Q free = 0, which automatically ensures that and by using Gauss's theorem (12) for the E field. In this way, D i V is expressed in terms of E and P averages: The average E i V can, in turn, be expressed in terms of E flux through a proper face of the unit cell. For instance, for the z-component one has is E flux through the cross section S z (z) of the unit cell defined by the plane z = const, see where Q(z) is the total induced charge within the volume V (z) enclosed by the boundaries of the unit cell and the plane z = const (see figure 4). Equation (33) is a consequence of Gauss's theorem (12) applied to the volume V (z) and accounting for the field periodicity. Now, accounting for (33) as well as E z (0) = const, V = a 3 , and that E z (0)/a 2 is nothing but E z S bott , we can write (32) as Similar equations for the x and y components can be obtained by replacing z → x or z → y so we may write S i being the face of the unit cell defined by the plane x i = 0. Note, as an insight, that (35) proves e. E field cannot be surface averaged in a flux-like manner. Now, let us show that the contribution to D i V due to the integral term appearing in (35) is completely compensated by P i V . From the definition of the average polarization where p total is the total dipole moment within the unit cell, we have The term dz S(z) ρ(r) dS is nothing but the total charge dQ(z) of the layer of thickness dz bounded by the planes z = const and z + dz = const (see figure 4). Therefore, Integrating the right-hand side by parts and accounting for Q(a) = 0 (recall that we consider lattices of initially uncharged inclusions) gives Similar expressions for the x and y components are obtained from this one by replacing z → x or z → y thereby giving By combining equations (31), (35) and (37), we finally obtain This result holds true for inclusions of arbitrary materials, either dielectric or metal. It follows from the above derivation that the surface S i appearing in (38) is the cell face defined by x i = 0. Due to the field periodicity, S i may be chosen on either of the two faces of the unit cell that are perpendicular to ith-axis. Unlike the derivation for the B field, we cannot, generally speaking, choose S i located anywhere between the faces. However, in case of dielectric inclusions one may still select S i arbitrarily since D flux along ith-axis is, like B flux, conserved in this case.

So far, we almost always treated the quantities D i V and B i V as if the local values D(r)
and B(r) are known. (The only exception was made in the preceding subsection.) However, electrodynamic solvers only usually provide information on local E and H fields, not D and B. Therefore, the question arises, how should D(r) and B(r) be calculated? One might calculate D and B from the known values of the local E and H fields, respectively, by using either multiplicative or additive relations [34] or where ε ≡ ε(r) = ε 0 ε r (r), µ ≡ µ(r) = µ 0 µ r (r); ε r (r) and µ r (r) are local values of the relative permittivity and permeability at point r; P = P(r) and M = M(r) are local polarization and magnetization. Surprisingly, despite being often considered equivalent, equations (39) and (40) predict different results for ε eff and µ eff , the results from (39) may even be paradoxical. To illustrate this, let us consider effective permittivity of a lattice of metal inclusions. From (39) it follows that Splitting the volume integrals into two parts corresponding to the region of the inclusion of a unit cell, V incl , and the region of the cell exterior to the inclusion, V ext = V − V incl , gives where it was accounted that ε r = 1 in V ext (vacuum or air outside the inclusions). For a metal inclusion in a static field, E = 0 inside the inclusion. Therefore, Then, the diagonal components of the effective permittivity of the lattice should be The obviously incorrect result ε eff = ε 0 has already been mentioned (yet for a particular case of thin metal inclusions) in [6]. In contrast, the additive relations (40) which, in fact, are definitions of D and B, lead to correct results. Direct averaging of (40) over the volume of a unit cell yields the average polarization and magnetization are related to the total electric and magnetic dipole moments of the unit cell, p total and m total : The moments p total and m total are due to charges and currents induced in the inclusion of a unit cell. It is these charges and currents that are responsible for the deviation of ε eff , µ eff from the corresponding free space values ε 0 , µ 0 . To summarize, the multiplicative constitutive relations, equations (39), should not be used in calculation of ε eff and µ eff of lattices of metal particles since they do not properly account for the contributions to D and B due to the charges and currents induced in a unit cell. 2. The two averaging procedures-the line averaging for (E, H ) fields and surface averaging for (D, B)-stem from different properties of these fields: the properties of E and H fields are governed by contour integrals (10) and (11), while those of D and B fields are governed by surface integrals (12) and (13). Note that the averaging formula for the D field is a consequence of Gauss's theorem for the E field rather than D itself. The reason for this is explained in section 2.2.4. 3. Periodicity of the total fields plays a crucial role in validity of the averaging formulae.
Maxwell's equations (10)-(13) themselves cannot guarantee that (volumetric) averages are equivalent to either line or surface averages. It is the field periodicity (14) which makes the equivalence possible.

Magnetic response of lattices of non-magnetic conducting inclusions (such as made of
Cu, Al, Ag, Au, GaAs, etc), being due solely to induced conduction currents, is only possible in a time-varying magnetic field that generates the currents by Faraday's law. In the purely static case we considered so far, the response is zero. Therefore, the approach based on Maxwell's integral equations for static fields cannot capture the artificial magnetism of such lattices. For an extension of this approach to time-varying fields, see the next subsection.

Modified field averaging
For a lattice in external time-varying electromagnetic fields the results of section 2.2 become, strictly speaking, incorrect. The reason for this is that the spatiotemporal variation of external fields causes the total fields in a lattice to be non-periodic and non-conservative. Accordingly, the conditions break down which are necessary for the replacements F i V → F i L and F i V → F i S to be valid. We may expect that the PHRS averaging formulae are still approximately valid in a quasistatic regime (say, at a < 0.1λ 0 ) when field periodicity holds approximately true and the circulation of total fields around the edges of unit cell faces is small. Beyond the quasistatic regime these formulae, although being previously claimed to work well [16], lead to implausible results (for numerical evidences, see our previous work [18] and section 3 below).
In order to extend the range of validity of the PHRS averaging formulae to appreciable unit cell sizes a ≈ 0.4λ 0 , a correction has been proposed [17] to account for the phase variation of the total fields across the unit cell. Alternatively, one could try to derive rigorous averaging formulae for non-static case by solving the respective boundary-value problem for dynamic fields in a unit cell. In the present paper, we obtain new averaging formulae somewhat heuristically: we start from the above derived formulae for static fields and then modify them to account for the phase variation of the fields across the unit cell.  (45). The average value E i of ith component of E field is calculated over any cross section S i of the unit cell that is parallel to the coordinate planes containing ith-axis.

Alternative averaging of static E and H fields.
Let us first note that in calculating averages of E and H fields, we actually are not bound to integration over line segments, see (6) and (7). In fact, integration over any surface within the unit cell that is parallel to the respective axis will give the same result. To show this, let us consider the volumetric average of z component of E field In section 2.2.1, we calculated this average as However, this choice of calculation is not unique: by using the constancy (18) of the line integral of E field, we can combine the three line integrals in (42) in two additional ways: either or The surfaces S x and S y do not necessarily coincide with the corresponding faces of the unit cell and form a set of planes parallel to the z-axis, see figure 5. Thus, they will be labeled as S z . The prime here is used to prevent confusing S z with the faces S z whose normal is parallel to the z-axis. The latter appear in the PHRS averages of D and B fields. Now, both equations (43) and (44) can uniformly be written as and for arbitrary component of E field we have Similarly one can obtain the respective formula for the H field average Hence, the averages of static E and H fields in a lattice can be calculated in three different ways-either as volumetric, surface or line averages, The value of equations (46) and (47) is in that they provide an alternative recipe for calculation of the average E and H fields. In the purely static case, it makes no difference: the equivalence · · · V = · · · S i = · · · L i of the three averages holds strictly true. Beyond the static case, one of the averaging methods may be more advantageous.

2.3.2.
Beyond statics: optimal field averaging and modified definitions of ε eff and µ eff . In a non-static regime, the phase advance of an incident wave leads to a spatial variation of its electric and magnetic fields across unit cells, see figure 6(a). While the surface averaging introduced by PHRS for D and B fields can account (within certain accuracy) for the phase variation of these fields, the line averaging introduced for E and H fields cannot. The latter uses only information on the local field values at the points with the same coordinate along the wave propagation directionk. Therefore, the integrals over two different lines-e.g. L 1 and L 2 in figure 6(b)-will give different average field values, even for empty space. In contrast, integration over the surfaces S i that are parallel to k (see S 1 , S 2 and S 3 in figure 6(b)) does account for the field variation. In the case of an incident plane wave, such S i are not necessarily the cell's faces: due to constancy of the plane wave fields in the directions ⊥k (see figure 6(a)), averaging over any cross section of the cell that is parallel to the faces (e.g. S 2 in figure 6(b)) returns the same result.
Based on this observation, the optimal averaging of E and H fields should employ integration of the local field values over surfaces that are parallel to the ith-axis and the wave vector k. Such surfaces form a subset of the previously introduced surfaces S i appearing in (46)-(48) and shown in figure 5. To differ these specific surfaces from S i , we add a subscript k to S i thus denoting them as S i,k . In this way, formulae of optimal averaging of E and H fields take the form and the formulae for effective parameters now read Expressions (50) H ) and (D, B) fields: the averages of the ith component of latter fields are computed over a plane that is normal to the ith-axis (flux-like surface averages); rather, the averages of E and H components are calculated over a plane that is along the ith-axis (in-plane or tangential averages). Formulae (50) for ε eff , µ eff together with the field averaging formulae (27), (38) and (49) are the basis for our method of homogenization of periodic materials. The advantages and the accuracy of this method are discussed below.

Numerical results and discussion
In order to illustrate our theoretical findings discussed in the preceding section 2, we have calculated the effective parameters of several simple SC lattices based on our equations (49) and (50). The calculated values were then compared with those obtained by the PHRS method. The difference between the two methods is expected to become substantial in a non-static regime, as is described in section 2.3.2. Also, the results from the two methods are compared, whenever possible, to the effective parameter values calculated by the volumetric averaging and TLEM.
The local field data used in the calculations were obtained from CST Microwave Studio Suite by simulating the propagation of a plane wave incident normally onto a slab of a periodic material. The total thickness of the slab was nine cells of size a each. The necessary field averaging was performed on the unit cell at the center of the slab.
The time-domain solver was used to calculate E and H field distributions inside the center cell and was set to terminate after the energy of the system had run down to -80 dB. Automatic meshing of the geometry was set to 50 lines per wavelength (solution convergence of this was verified but is not presented here). Once calculated, post-processing templates within the software were used to calculate the necessary line, surface and volumetric averages of E and H fields. Since CST solvers do not solve for either (D,B) or (P,M) fields, local values of D and B fields were calculated from those of E and H using the multiplicative constitutive relations (39). Averaging of (D,B) fields was then done following the same post-processing as for (E,H ) fields. Ratios of the averaged quantities were finally used to calculate the effective material parameters as functions of the unit cell size. In the case of conducting spheres, the use of multiplicative relations (39) renders calculation of volumetric averages incorrect (see section 2.2.5) so results from volumetric averaging were not calculated for this case and are not provided below.
The examined region of the unit cell sizes extends up to a = (0.25 − 0.35)λ 0 . At larger sizes, the resonance and diffraction effects become important and the very concept of metamaterials homogenization is expected to fail.

Empty cells
The case of free space not only provides the simplest possible verification of the theoretical results obtained in previous sections but also gives an insight into how field averaging works for non-empty lattices. Additionally, this case is very instructive since some conclusions can be (and, actually, have been) drawn as to possible corrections one has to make to the initial form of the field averaging originally proposed by PHRS.
The free space permittivity and permeability can be treated as effective parameters of a lattice of empty cells and hence can be calculated by using either the PHRS or our averaging formulae. The calculation results are presented in figure 7.
As seen from figure 7, the effective parameters ε eff , µ eff of empty cells calculated by the PHRS method (solid black curves) are complex valued and behave in an oscillating manner. Their real and imaginary parts start from correct values ε eff,r = µ eff,r = 1, ε eff,r = µ eff,r = 0 at a/λ 0 = 0 and then deviate more and more as the unit cell size increases to a = 0.75λ 0 (for the real parts ε eff , µ eff ) and a = 0.5λ 0 (for the imaginary parts ε eff , µ eff ) and asymptotically approach zero in the limit a → ∞. In contrast, our method (dashed blue curves) gives exactly ε eff,r = µ eff,r = 1, ε eff,r = µ eff,r = 0 regardless of the unit cell size.
The oscillating behavior of the effective parameters of empty cells within the PHRS approach was first discussed by Smith et al [17] where the following analytical expression was derived for µ eff : (k is the wave number of a plane wave E 0 exp[i(k · r − ωt)] propagating along one of the axes of the unit cell). A similar expression can also be obtained for ε eff , so the relative effective permeability and permittivity of empty lattice are Equation (52) suggests that both permittivity and permeability of free space are spatially dispersive and complex valued, their real and imaginary parts being It is these results that are accurately reproduced by the solid curves in figure 7.
Later, Smith and Pendry [16] obtained a real-valued expression for ε eff and µ eff , which contradicts (52) but is similar to its real part (53). (Note that Smith and Pendry labeled the lattice period 2a so in their notations the factors of 2 do not appear in (55) making it of completely same form as (53).) Expression (55) may be considered as a somewhat better alternative to (52) as it predicts ε eff = µ eff = 0, although still predicting the oscillating behavior of ε eff and µ eff . However, a thorough examination reveals that 1. both the expressions (52) and (55) are particular cases of more general ones that still retain both inconsistencies-the oscillating behavior and complexity of ε eff and µ eff ; 2. the difference between (52) and (55) is merely due to different choices of the integration path L that have implicitly been made in [16,17] while deriving these expressions.
The more general expressions for ε eff , µ eff can be obtained, within the PHRS method, by calculating the averages E i L , H i L for an arbitrary location of the appropriate path L and read Here it is assumed that L is in the direction of E field of the incident wave (see (20)) and, thus, L⊥k (see figure 6), and l is the coordinate of L in the k direction (0 l a). Results (52) and (55) are, then, particular cases of (56) corresponding to l = 0 and a/2. Note that according to (56), ε eff , µ eff are not only complex, oscillating and spatially dispersive, but also dependent on the location of the path L. The latter is illustrated in figure 8.
As seen from figure 8, the effective parameters of free space calculated by the PHRS method become real only when the integration path is located in the middle of the unit cell face. The choice l = a/2 can thus be considered optimal for this method, although it does not comply with the PHRS prescription to use l = 0 [6]. Note, however, that even with this optimal choice, the free space parameters are still oscillating, see (55) and the corresponding dashed red curve in figure 8(a).
Physically, these issues with ε eff , µ eff are caused by the spatial variation of the electric and magnetic fields of the incident wave over the empty cells. Formally, they are due to the oscillating sinc and sinc-like factors [sinc(x) ≡ (sin x)/x] appearing in equations (52)-(56). The factors stem from the line integrals defining the E and H averages within the PHRS approach.
To avoid the issues and match the values of ε eff , µ eff calculated by the PHRS method with the true values ε 0 , µ 0 of free space, it has been proposed [16,17] to simply remove the oscillating factors from the expressions for ε eff , µ eff . Such 'handmade correction' was also suggested to give plausible results when applied to non-empty lattices. We shall examine this correction in the subsections that follow.
In this subsection, we have computed effective parameters of SC lattices of homogeneous spherical inclusions. Dielectric spheres of two different kinds were considered-non-magnetic and magnetic. In contrast to Smith and Pendry [16], who considered only lattices with small fill factors f 0.065, below we provide data for f = 0.2, 0.3.
Seven different methods were employed in our calculations: four versions of the PHRS method, our proposed method, TLEM and conventional volumetric averaging. We specify which version of the PHRS method is being indicated by adding additional labels such as l = 0 (integration line L chosen at the edge of the unit cell), l = a/2 (integration line L halfway along the face of a cell) or 'corrected'. The correction, whenever applied, was performed by dividing computed complex values of ε eff , µ eff by the factor equal to the right-hand side of (56), with k = k stat ≡ 2πn stat eff /λ 0 = (ε stat eff,r ) 1/2 ω/c, where ε stat eff,r is the relative effective permittivity calculated in the static limit (i.e. at a → 0 or, equivalently, ω → 0).

Non-magnetic spheres.
We considered a lattice of spheres made of a non-magnetic material with ε r = 8, µ r = 1. Figure 9 shows our numerical results for both the real and imaginary parts of the effective permittivity of the lattice. The choice of l value and whether the correction has been applied are both indicated in the legend. Since non-magnetic dielectric spheres possess a weak magnetic response, below we do not analyze their permeability.
As seen from the figure, the results from all the seven methods converge to each other in the long-wavelength limit. This confirms the conclusion we made in section 2 on equivalence, in the static regime, of the PHRS method, volumetric averaging and our method. Conversely, as the unit cell size increases, the results diverge. Note, first of all, that both the initial PHRS method and its corrected version (see, respectively, the black dashed and solid curves in panels To quantify the difference between the methods we consider, we have calculated a relative deviation ε eff ≡ |ε eff − ε vol eff |/ε vol eff of the values ε eff calculated by different methods from the value ε vol eff calculated by volumetric averaging. The deviation versus unit cell size is shown in panel (c) for the most accurate methods presented in the inset. As seen from figure 9(c), at a/λ 0 = 0.35 our method provides a substantially more accurate value of the effective permittivity when compared to the optimized corrected version of PHRS method ( ε eff = 6% versus 12%). Also, up to a/λ 0 ≈ 0.25 our proposed method returns effective permittivity values that are nearly identical (with ε eff < 0.5%) to those obtained by volumetric averaging.
The results presented in figure 9 suggest that PHRS method is highly sensitive to the choice of the integration path location not only in the trivial case of empty cells (see the preceding subsection) but also in the case of non-empty lattices. Moreover, the initial (with l = 0) version of this method and its corrected modification both provide, similar to the case of empty cells, complex valued effective permittivity, despite the fact that the inclusions we consider here are lossless.
Furthermore, it turns out (the results are not shown in this paper) that the first three of the above considered modifications of PHRS method predict implausible behavior of both ε eff and µ eff for inclusions of various sizes, shapes and materials. Based on this, we remove these modifications from further discussion and shall keep only the corrected PHRS method which employs the optimal choice l = a/2.
Finally, it is worth noting that, according to figure 9, the effective permittivity of dielectric spheres disperses monotonically upward in the region a/λ 0 < 0.35. In terms of frequency dependence of the effective index of refraction n eff (ω) = [ε eff,r (ω)µ eff,r (ω)] 1/2 this means (accounting for µ eff ≈ 1) that the lattice has normal dispersion at ω < 0.35 × (2πc/a) ≈ 2.2c/a. Normal dispersion of ε eff , µ eff and n eff of lossless lattices beyond resonance regime is well known (see e.g. [39]) and has been confirmed numerically by different methods (e.g. in [37,40]). In contrast, Alù [11] calculated effective parameters of dielectric spheres within his homogenization theory and obtained a different behavior of ε eff with anomalous dispersion at low frequencies, see figure 3 in [11].

Magnetic spheres.
Lattices of particles made of magnetic non-conductive materials such as ferrites, magnetic polymers, magnetodielectric composites and spin glasses may have a substantial magnetic response. Under proper tuning of their effective permittivity and permeability, such lattices can be used in the design of left-handed metamaterials [35,41]. In figure 10, calculation results are presented for effective permeability of an SC lattice of spheres having ε r = 3, µ r = 4, at lattice fill factor f = 0.2. The results were obtained by using the four most accurate methods mentioned in the preceding subsection.
One can see from figure 10(a) that all the four methods predict the same static value of the effective permittivity. As the unit cell size increases, our proposed method, TLEM and volumetric averaging give values of µ eff,r that are qualitatively close to each other. In contrast, the value calculated with the optimal PHRS method diverges substantially from these three. For example, at a/λ 0 = 0.3 the relative deviation µ eff ≡ |µ eff − µ vol eff |/µ vol eff is less than 2% for our method and about 12% for the PHRS method, see figure 10(b).
Based on figures 9 and 10 we may conclude that for lattices of non-conducting particles our proposed method provides better results than any modification of the PHRS method does.

Perfectly electrically conducting (PEC) spheres
Lattices of metal spheres have long been considered both theoretically and from the viewpoint of applications dating back to, at least, the works by Rayleigh [1] and Kock [42], respectively. In general, lattices of conducting inclusions may possess appreciable both electric and magnetic  response since the individual particles can easily gain substantial polarization and magnetization in time-varying external fields. Therefore, the lattice permittivity and permeability are both of interest in this case. Figure 11 shows calculation results for the real parts of the effective parameters of an SC lattice of perfectly electrically conducting (PEC) non-magnetic spheres at lattice fill factor f = 0.2. Volumetric averaging results are not shown here for reasons discussed in the beginning of section 3. Therefore, the quantities ε eff and µ eff shown in figures 11(c) and (d) present deviations from the corresponding TLEM values. The data are calculated for unit cell sizes a 0.25λ 0 . We were unable to extract the effective parameters with TLEM at larger sizes because of strong electromagnetic wave scattering observed at a > 0.25λ 0 .
As seen from figures 11(a) and (b), the results from all the methods, again, converge in the long-wavelength limit to the same static values of ε eff and µ eff . The static permittivity value ε eff,r = 1.756 calculated here by three numerical methods has been obtained previously using a different numerical technique [43] and from analytical theories [4,40].
Beyond the static regime the results obtained by PHRS and our method diverge away from those calculated via TLEM. A satisfactory agreement between the three methods exists only at smaller unit cell sizes. For instance, deviations ε eff , µ eff of less than 5% can be achieved at a 0.15λ 0 , as seen from figures 11(c) and (d). Therefore, one may conclude that for conducting inclusions averaging over the boundaries of unit cells gives less accurate results as compared to dielectric inclusions.
Note, again, that all the three methods presented in figure 11 predict an upward dispersion of ε eff and a downward dispersion of µ eff . The same behavior was reported, e.g. in [44,45] while Alù [11] obtained an opposite one with a negative slope for ε eff and positive for µ eff , see his figure 9.

Conclusions
Homogenization of metamaterials by field averaging over the boundaries-edges and faces-of their unit cells was proposed by PHRS as an alternative to conventional volumetric averaging introduced by Lorentz. The latter was believed to fail in case of lattices of thin metal inclusions such as wires and SRRs [6,16]. In contrast, PHRS method was suggested to be valid for lattices with virtually no restrictions on the contents and size of the unit cell [16]. We have shown that PHRS averaging formulae can, in fact, be rigorously derived from Lorentz definition of the field averages. However, volumetric field averages F i V can be reduced to corresponding PHRS averages F i S or F i L only in the case of static fields. For time-varying fields such a reduction is, strictly speaking, impossible. Therefore, the PHRS homogenization scheme can be considered merely a special case (approximately valid in quasistatic regimes) of the more general Lorentz scheme.
Beyond static field, the PHRS method not only predicts complex valued effective parameters even for lossless lattices but also becomes sensitive to the choice of the path L of integration of curl-conforming fields E and H . In this case, plausible values of the effective permittivity and permeability can be obtained only under specific choice of L and, additionally, if the correction due to Smith et al [17] is applied. As a special note, the correction factor approach proposed by Smith et al [17] also assumes a non-dispersive effective material. As dispersion increases, the accuracy of the correction factor decreases.
Alternatively, one can obtain plausible results by employing surface averaging of E and H fields instead of the line averaging suggested by PHRS. We have theoretically proven that averaging of these fields over properly chosen faces of the unit cell is fully equivalent to either volumetric or line averaging in the static regime. Although the equivalence deteriorates in time-varying fields, the surface averaging can still account for the phase variation of the total fields across the unit cell and thus provides better results automatically without invoking any correction procedures.
Overall, homogenization of periodic metamaterials by field averaging over unit cell boundaries can be considered a good alternative to Lorentz averaging or TLEM in quasistatic regimes (a < 0.1λ 0 ). Particularly, it can provide highly accurate static values of ε eff , µ eff regardless the geometry and material of the inclusions. At the same time, it uses only information on the local fields E(r) and H(r) on the boundaries of the unit cell thus substantially decreasing the required computational resources and enabling computation of the effective parameters even if information on P(r) and M(r) is not available.
For larger unit cell sizes, averaging over cell boundaries should be performed with care. In intermediate regimes-at a = (0.1 − 0.35)λ 0 -the accuracy of 2-20% for typical cases shown in the present paper can be achieved, less accurate results being expected for lattices of inclusions with high dielectric contrast. At a 0.4λ 0 , the resonance and diffraction effects become crucial so all field averaging methods are expected to fail as well as the very concept of metamaterial homogenization. Although in this regime one can still formally calculate the effective parameters by using a variety of field averaging techniques and the local field data provided by a solver, the homogenized medium cannot describe the properties of the highly scattering simulated metamaterial.