Boundary conditions for quadrupolar metamaterials

One of the long-standing problems in effective medium theories is using the knowledge of the bulk material response to predict the behavior of the electromagnetic fields at the material boundaries. Here, using a first principles approach, we derive the boundary conditions satisfied by the macroscopic fields at interfaces between reciprocal metamaterials with a quadrupolar-type response. Our analysis reveals that in addition to the usual Maxwellian-type boundary conditions for the tangential fields, in general—to ensure the conservation of the power flow and Lorentz reciprocity—it is necessary to enforce an additional boundary condition (ABC) at an interface between a quadrupolar material and a standard dielectric. It is shown that the ABC is related to the emergence of an additional wave in the bulk quadrupolar medium.


Introduction
Effective medium theories are extremely useful to describe 'low energy' wave phenomena, as they typically allow the reduction of complex propagation problems in heterogeneous media to the wave propagation in 'continuous media' described by a few effective parameters. For example, the interaction of electromagnetic waves with natural media is typically described by a permittivity and a permeability, which depend on the electronic properties of the atoms [1].
Similarly, in semiconductor physics the electron transport properties may be described by an effective mass [2].
While effective medium methods are quite common in many branches of physics, the topic of determining the effective medium response and of using that knowledge to solve relevant physical problems remains to this day challenging and is characterized by several unsolved and unclear issues. The difficulties are well illustrated by the problem of wave propagation in electromagnetic periodic media (metamaterials).
The first difficulty is concerned with the very definition of the effective parameters for bulk (unbounded) structures. Evidently, effective medium approaches provide necessarily an approximate description of the wave phenomena, and thus depending on the assumptions that are made one can obtain different effective parameters. Thus, one can find in the literature different methods for calculating the effective response, which rarely yield the same results, even in the long wavelength limit [3][4][5][6][7][8][9][10][11][12][13]. Throughout this work, we adopt the effective medium framework originally introduced in [6, [14][15][16][17], based on the idea of source-driven homogenization.
The second difficulty-which is the research topic of this article-is more subtle, and concerns using the bulk effective response to describe the wave propagation in the presence of interfaces between different media. The trouble is that the macroscopic electromagnetic fields can vary abruptly at material boundaries, and characterizing the jump discontinuities of the fields can be a quite formidable task. This theme is of obvious and fundamental importance because the most interesting wave phenomena occur at the interfaces between different materials. Solving the problem in the most general case seems to be a hopeless task, particularly if the medium response has a strong spatial dispersion. Indeed, the fields' behavior near the boundaries typically depends on the particulars of the microstructure, and thus in general the knowledge of the bulk responses may be insufficient to characterize the interface response.
The goal of this article is to derive the boundary conditions for the macroscopic electromagnetic fields at interfaces of media with a weak spatial dispersion, in the 'quadrupolar media approximation'. The quadrupolar media approximation is based on the assumption that the response of the 'meta-atoms' is equivalent to that of a collection of electric dipoles, magnetic dipoles, and electric quadrupoles. Here, using first principles arguments it is shown that this assumption about the microstructure is sufficient to determine some boundary conditions satisfied by the electric and induction fields at material interfaces. There is an important body of work on the topic of boundary conditions for materials with a quadrupolar response, pioneered by Raab and co-authors, with more recent contributions by Yaghjian [18][19][20][21]. The analysis of these authors is based on the macroscopic (homogenized) Maxwell's equations. Interestingly, even though our theory is in essence microscopic, it is shown that it gives results compatible with previous studies [18][19][20][21]. However, for the first time we highlight the need for an additional (non-Maxwellian) boundary condition (ABC), and link the ABC with the emergence of additional waves in the bulk medium. We prove that in the simplest case (if the macroscopic electromagnetic fields do not have δ-type singularities), the ABC at the interface between a quadrupolar material and a standard dielectric isˆ⋅ ⋅ = Q n n 0,n being the unit vector normal to the interface and Q the quadrupole moment density. Different aspects of the quadrupolar response of metamaterials and boundary conditions for media with weak spatial dispersion have been discussed in other works (e.g. [22][23][24][25]), but the topic is clearly in its infancy and is largely virgin territory. This article is organized as follows. In section 2 we set the stage for the analysis that follows and briefly review the electrodynamics and homogenization of quadrupolar media. In section 3 we derive the boundary conditions for interfaces of reciprocal quadrupolar media using both microscopic and macroscopic approaches. It is shown that the macroscopic fields must satisfy an ABC at the interface of a quadrupolar material and a standard dielectric. In section 4 we outline how the effective parameters associated with the quadrupole response can be determined, and the ABC is linked to the emergence of additional waves in the bulk quadrupolar medium. Next, in section 5 a numerical example is used to illustrate the application of the developed methods to a metamaterial formed by pairs of metallic nanorods. A summary of this study is given in section 6.

Quadrupolar media
We consider the problem of wave propagation in a generic electromagnetic metamaterial. The macroscopic Maxwell's equations are obtained after averaging the 'microscopic fields', and for metal-dielectric metamaterials can be written as [1,26]: where E, B and j e , are the macroscopic electric field, the macroscopic magnetic induction field, and the macroscopic external current density, respectively. The generalized polarization vector is determined by the spatially averaged microscopic current density , where ε ε = r ( ) denotes the permittivity of the inclusions and e the microscopic electric field [6,26]. We follow Russakoff [1,27] and suppose that the macroscopic fields are obtained after weighting the microscopic counterparts with a suitable test function f. Such an averaging procedure preserves the structure of Maxwell's equations. In particular the generalized polarization vector P g is given by: As is well known, the convolution operation is equivalent to spatial-filtering. Following our previous works [6,17,26] otherwise. This type of averaging corresponds to ideal low-pass spatial filtering of spatial harmonics in the BZ.
To obtain a multipole expansion for P g [1,28], we suppose that the microstructure of the material is such that each 'inclusion' can be well approximated by a collection of multipoles centered at the inclusion's geometrical center. This is equivalent to stating that:  I  I  I  I  I  I  c   d , ,

,
In the above, r I is the position vector of the Ith particle in the material, ∂ = ∂ ∂x / n n , and the vectors P I , Q n I, , R nm I, determine the multipole moments of the Ith particle normalized to the volume V c of the unit cell. Because the operation of averaging commutes with the spatial derivatives [1,27], we have: The symbol '*' represents the operation of convolution, and f is the test function. We want to obtain explicit formulas for the multipole moments in terms of the microscopic current. P I is obtained by integrating the microscopic current (3) over the volume ( )  and higher multipole moments are negligible. This is the 'quadrupolar medium approximation'. Substituting the previous formulas into equation (5), one obtains the following explicit definitions for P and= ∑ ⊗ Q u Q n n n (the symbol ⊗ represents the tensor product of two vectors): Evidently, the vector P is the usual polarization vector. The tensor Q is designated here as the total quadrupole density. As outlined in appendix A, it includes the averaged response of both the electric quadrupoles and of the magnetic dipoles, such that el where I is the identity dyadic, Q el is the standard electric quadrupole moment density and M is the standard magnetization vector [1]. Note that the total quadrupole density determines univocally the electric quadrupole density and the magnetization vector. Specifically, Q el is the symmetric part of the tensor Q , whereas M is determined by the anti-symmetric part. The relation between Q el and M and the microscopic currents is reviewed in appendix A. Within the quadrupolar approximation, the multipole expansion (4) becomes: The divergence of a dyadic is by definition , whereû i is a unit vector along a coordinate axis.
It is well known that a decomposition of the type  = − ⋅ Q P P g is not unique. However, there is no arbitrariness in our formulas since they follow naturally from the hypothesis that to a good approximation the response of the microscopic current can be approximated by point multipoles (equation (3)). Therefore, our P and Q have strict physical meaning attached to them. Obviously, for quadrupolar media, the Maxwell's equations reduce to: For future reference, we note that in the case where the microscopic electromagnetic fields have the Bloch property, with propagation factor ⋅ e k r i being ∈ B Z k . ., the wave vector, the polarization vector reduces to: In the above, r 0 is the center of the inclusion in the unit cell Ω ( ), and G is a generic vector in the reciprocal lattice. For simplicity, it was assumed that the unit cell contains a single inclusion (alternatively, all the inclusions may be regarded as a single one). Similarly, for Bloch waves the quadrupole density and the generalized polarization satisfy, ⎛ r k r r Thus, P, Q , and P g are plane waves, and their amplitudes are obtained by suitable integration of the microscopic currents over the unit cell.
In the next section, we highlight the relevance of the effective response of P and Q to an external excitation, and make clear that the material response associated with Q is essential to characterize the electrodynamics of the macroscopic fields in problems that involve interfaces of quadrupolar media.

Boundary conditions from a microscopic approach
Our objective is to characterize the boundary conditions satisfied by the macroscopic electromagnetic fields at an interface between quadrupolar media. Previous works have analyzed this problem based on macroscopic approaches, which essentially are variants of the standard technique of enclosing the interface with a 'pillbox' of infinitesimal volume and integrating the macroscopic Maxwell's equations [20,29,30]. Here, we use instead a microscopic theory adapted from our previous work [31]. The findings of this subsection provide an alternative derivation of results previously reported in the literature [18][19][20][21], and particularly they confirm the boundary conditions originally discovered by Raab and co-authors [18, p 156].
In [31] the concept of transverse averaged (TA) fields was introduced to study in a systematic manner the scattering by metamaterial screens. The TA fields are defined as the microscopic fields averaged only along the directions parallel to the pertinent (planar) interface. Thus, a generic TA field is given by the convolution of the associated microscopic field with a test function f TA , similar to the approach followed to define the bulk fields [1,26], but with the important nuance that f TA is of the form , δ being the Dirac distribution. For example, for an interface parallel to the xoy plane the TA electric field is given by: TA TA where e represents the microscopic electric field. Note that the microscopic field is averaged exclusively along the transverse coordinates. As in section 2, it is assumed that f TA is such that the convolution operation is equivalent to an ideal low-pass filtering along the x and y directions. In case of a periodic structure, the 'bandwidth' of the ideal low pass filter is determined by the transverse Brillouin zone. In contrast, the usual bulk fields E and B used in effective medium theories are averaged over a volumetric region, i.e. they are also averaged along the z direction normal to the pertinent interface. It is straightforward to derive the boundary conditions satisfied by the TA fields at an interface [31]. For example, the tangential components of the TA electric and induction fields, E TA and B TA , are continuous at the interface (there are some exceptions discussed in [31] but these are irrelevant to this work). It was shown in [31] that for media which at the microscopic level consist of a collection of electric and magnetic dipoles, it is possible in the long wavelength limit to write E TA and B TA in terms of the bulk fields E and B, and in this manner find the boundary conditions for the bulk fields at the interface. Here, we generalize these ideas to the case of quadrupolar media. It is important to highlight that the boundary conditions for E TA and B TA are different from the boundary conditions for E and B, and in this sense one can say that the boundary conditions depend on the considered effective medium theory.
For conciseness, the technical details are given in appendix B. Our analysis shows that at an interface between two generic quadrupolar media the tangential macroscopic fields must satisfy: wheren is the unit vector normal to the interface, Grad is the surface gradient, and= ⋅ ⋅ Q Q n n nn , and ⌊ ⌋ = − + − F F F represents the jump discontinuity of a generic vector (or dyadic) at the interface.
The boundary condition (13b) is a generalization of the usual Maxwellian boundary condition for the magnetic field. In fact, when the electric quadrupole vanishes, which is tantamount to stating that the tangential component of is continuous at the interface, consistent with classical theory.
On the other hand, equation (13a) generalizes the Maxwellian boundary condition . Note that somewhat surprisingly, (13a) allows for a macroscopic electric field discontinuous at the boundary. This will be discussed further below.
It is straightforward to obtain the jump conditions for the normal components of the electromagnetic field. Applying the surface divergence operator (Div) to equations (13b) and (13b) and using the Maxwell equation (10) it is readily found that, 0 g tan where  = − ⋅ Q P P g and the subscript 'tan' denotes the projection of a vector onto the interface. The jump conditions for the normal fields are in some sense redundant because they follow directly from the boundary conditions for the tangential fields and the Maxwell's equations. When Interestingly, the boundary conditions (13) derived with our microscopic approach are completely consistent with what was obtained in [18, p 156] (see also [20,21]) and [19] using a macroscopic theory. Thus, up to now our findings basically confirm earlier work.

An additional boundary condition for sectionally continuous macroscopic fields
Next, we consider again an interface between two materials, and suppose that the macroscopic electromagnetic fields, the polarization vector, and the total quadrupole density are sectionally continuous in all space, i.e. δ-Dirac type singularities at the interfaces are not allowed. At first sight, this assumption may look obvious and one may be tempted to reject the possibility of having field singularities at the interfaces. However, in the vicinity of the interfaces there is a thin transition layer where the bulk constitutive relations usually break down [19]. The effect of the transition layer can be modeled with equivalent sources placed at the interface, and the presence of these sources may cause the electromagnetic fields to be singular at the boundary [19]. Thus, it is unclear if δ-type singularities may be physically possible in quadrupolar metamaterials. In section 3.4 we further discuss this exotic possibility.
Under the assumption that the fields are sectionally continuous, the boundary conditions at interfaces of quadrupolar media can be readily obtained using the standard macroscopic approach of integrating the fields over a pillbox with infinitesimal volume. Indeed, integrating both members of the Maxwell equation (10) over a volume V it is readily found that: Thus, if V has infinitesimal volume and Q E B P , , , are sectionally continuous (δ -function-type singularities are not allowed) then: ⎛ Taking V as a pillbox surrounding the interface of interest, it is seen that in order that the surface integrals vanish it is necessary that: The second boundary condition in (18) is equivalent to the following two: After some algebra, it can be checked that the former condition yields equation (13b), while the latter is ⎢ ⎣ ⎥ ⎦ = Q 0 nn . On the other hand, the first boundary condition in (18),ˆ× ⌊ ⌋ = n E 0, only requires that E B , are sectionally continuous at the interface, and implies that in such a case the tangential component of the electric field is necessarily continuous. This is compatible with equation In particular, if δ-type singularities in E B , are rejected as assumed here, it is necessary that = Q 0 nn at an interface between a standard dielectric and a quadrupolar medium.
In summary, this analysis shows that if the macroscopic fields are sectionally continuous, the pertinent boundary conditions at the interface between a quadrupolar material and a standard dielectric are equation (13b) and nn Equation (20) may be regarded as an additional boundary condition (ABC) [32][33][34][35][36][37][38][39], which of course is not found in the classical Maxwellian theory. To our best knowledge, ABCs for quadrupolar media were not discussed in previous works (e.g. [18,19] is not really surprising, because spatially dispersive media are characterized by internal degrees of freedom whose dynamics are not described by the traditional boundary conditions [32,35]. Moreover, media with a nontrivial Q el may support an 'additional wave' (e.g. isotropic spatially dispersive media may support a longitudinal wave besides the usual transverse waves [40]), and hence the usual boundary conditions for the tangential components of the fields are insufficient to characterize the electromagnetic response. This statement will be made precise in sections 4 and 5.
Several recent works demonstrated that ABCs are essential to determine the response of metamaterials with strong spatial dispersion [36][37][38][39]. Here, we propose the ABC (20) to characterize interfaces of media described by an electric quadrupole response, i.e. characterized by weak (second order) spatial dispersion. Evidently, (20) is not equivalent to Pekar's boundary condition, which imposes that the component of the polarization current normal to the interface vanishes [32]. It should be mentioned that a previous study claimed that ABCs are a 'historical mistake' and unnecessary [34]. However, such a conclusion relies on the assumption that the fields in a finite thickness spatially dispersive slab are coincident with the fields radiated by a particular localized current distribution placed in the bulk medium. Unfortunately, [34] does not take into account other excitation forms that are allowed when an interface is formed.
In the framework of this subsection, the boundary conditions for the normal field components can be obtained by applying Div to equations (13a) (with = Q 0 nn ) and (13b), and using the Maxwell equation (10). This procedure yields exactly the same result as in equation (14). It is instructive to note that when = j 0 e Maxwell's equations (10) imply that , rather than to the result in equation (14). The flaw in the reasoning is that the standard pillbox approach can be used only when is sectionally continuous near the interface [29]. In the current framework, Q E B P , , , are sectionally continuous, but  ⋅ Q can have δ-type singularities at the boundary. This invalidates the formula

Lorentz reciprocity
We are interested in composite media formed by dielectric or metallic inclusions at the microscopic level. Clearly such structures are reciprocal, and this poses important constraints on the effective medium response. The analysis of this subsection is not restricted to the case wherein the fields are sectionally continuous at the interfaces, and thus our only assumption is that the boundary conditions (13) (derived with the microscopic approach) hold. Next, it is shown that in order that these boundary conditions are compatible with the reciprocity theorem the material response needs to satisfy additional consistency conditions. To this end, we consider a quadrupolar material region surrounded by a regular dielectric with permittivity ε ext and permeability μ μ = 0 g where ′ ↔ ″ represents the same expression as in the left hand side of the equation but with the superscripts ′ and ″ interchanged. (9)) it is found after simple manipulations that, i i i g and hence equation (21) is equivalent to: The next step is to integrate both sides of the equation over the volume V of the material. Using the divergence theorem it follows that, ⎡ wheren is the unit vector normal to the boundary surface ∂V ( ) of the considered region (figure 1), and all the fields are evaluated at the inner side of the surface. In the above, we used the identity, . Because the volume enclosed by ∂V is formed (at the microscopic level) by reciprocal structures, the following equation must be satisfied in the dielectric region (exterior to V): ⎛ where the subscript 'ext' indicates that the fields are evaluated at the outer side of ∂V . Subtracting member by member the above expression and equation (24), and using the boundary conditions (13), we see that the electromagnetic fields must satisfy: where we used the shorthand notation ⎛ g 0 Integrating by parts the first term of the surface integral it is found after rearranging the terms that: This is the consistency condition that we looked for at the interface between a quadrupolar medium and a standard dielectric. Specifically, the response of a reciprocal quadrupolar medium body (figure 1) must be such that the above condition is satisfied by any electromagnetic field distributions originated by sources placed outside the body.
Interestingly, with the additional hypothesis that the fields are sectionally continuous at the interface, the surface integral vanishes because the ABC (20) holds. To ensure that the volume integral is also zero, one can impose that the respective integrand vanishes: i i i We will see in section 4.3 that the above equation leads to the familiar reciprocity restrictions on the bulk material parameters. In summary, the ABC (20) together with condition (29) ensure that the Lorentz reciprocity theorem is satisfied. More generally, without limiting ourselves to sectionally continuous macroscopic fields, we see that to comply with the Lorentz reciprocity theorem it is enough that both equation (29) and  (28)). In the next subsection, we study in detail the subclass of quadrupolar materials that satisfy equation (30). It is proven that the electromagnetic fields in such materials must obey a generalized ABC at an interface with a conventional dielectric.

A generalized ABC at an interface with a standard dielectric
Next, we consider a reciprocal quadrupolar material with electromagnetic response consistent with equations (29), (30). Equation (30) can be rearranged to give: nn 0 g 0 Since the two field distributions are arbitrary we can have an identity only if both members are equal to a constant independent of the excitation. Hence, we can write: where α Q is a scalar that may depend on the particulars of the interface and on the frequency, but not on the excitation. Equation (32) generalizes the ABC (20) to the class of materials with a response consistent with equations (29), (30). It should be evident that for materials in this class the macroscopic fields at the interfaces may be singular. Indeed, if the scalar α Q is different from zero then ≠ Q 0 nn at the interface, and from section 3.2 it follows that either E or B have a δtype singularity. Obviously, by setting α = 0 Q the generalized ABC reduces to equation (20). For standard Maxwellian media with = Q 0 el , the macroscopic fields are sectionally continuous and thus α = 0 Q for such materials. As discussed previously, even though it is quite tempting to assume that the fields cannot be singular and that α = 0 where the subscript 'ext' denotes the fields evaluated at the dielectric side of the interface. Thus the generalized ABC (32) can be replaced by: The constant α Q -if different from zero-must be determined with a microscopic approach and depends on the specific quadrupolar medium.
In appendix C, we demonstrate that the generalized ABC is consistent with the conservation of the power flow at an interface provided α =

{ }
Re 0 Q . Moreover, it is shown in appendix C that the complex Poynting vector is given by: c 0 It is crucial to note that the boundary conditions (13) and (33) are written in terms of Q . Thus, it is essential to know the individual response of Q to an excitation to characterize the dynamics of the macroscopic fields near the interfaces. Note that knowledge of  = − ⋅ Q P P g is insufficient to determine Q . Interestingly, even though the total quadrupole density Q is written in terms of Q el and M (equation (8)), the individual responses of the electric quadrupole density and of the magnetization are not required to formulate the boundary conditions. Hence, in this sense the quantities Q el and M are of secondary importance as compared to Q .

Weak (second order) spatial dispersion
So far, we have not discussed either how the effective response of P and Q can be determined or what the associated constitutive relations are. The objective of this section is specifically to outline how this can be done, and to highlight that the requirement for an ABC at an interface between a quadrupolar medium and a standard dielectric emerges naturally in the case where the effective medium is characterized by spatial dispersion of second order.

Effective response
The key and most challenging aspect in homogenization theory is to establish a relation between the spatially averaged multipole moments and the macroscopic fields. Next, we show that the electromagnetic response of the polarization vector P and of the total quadrupole density Q can be determined by a straightforward generalization of the ideas of our previous works [6,[14][15][16][17].
The homogenization approach developed in [6,16] permits determining the response of P g to the macroscopic field such that in the Fourier spatial domain: , n = x,y,z) to the electric field is generally of the form: It is useful to note that because of equation (

Constitutive parameters for weak spatial dispersion
The responses of P and Q to the macroscopic electric field may generally depend on the wave vector in quite a complicated manner. Here, we are interested in the case of weak spatial dispersion such that ε P and q n are polynomial functions of k with degree equal to one or less. This is equivalent to assume that the following Taylor expansions hold: Substituting these formulas into equation (37), one sees that the global response ε ω k ( , ) ef is a polynomial in k of degree equal or inferior to two: Thus, this approximation allows for spatial dispersion effects of second order.

Lorentz reciprocity again
Next, we investigate the reciprocity restrictions implied by equation (29) on the effective medium parameters. Putting  = − k i in equation (38) we obtain after some algebra, where ∂ = ∂ ∂x / l l . It is evident that in order that the symbols ′ and ″ can be interchanged for arbitrary fields it is required that (below it is implicit that = k 0): These are the reciprocity constraints on the effective medium parameters obtained from equation (29). It can be readily verified that these constraints imply that the effective dielectric function given by equation (39)   where the tensors ζ μ , EB may depend on frequency, and η μ ε = / 0 0 0 is the host impedance. The tensor μ represents the magnetic permeability and the dimensionless tensor ζ EB is associated with magneto-electric coupling effects as further discussed below. Under these hypotheses, it is straightforward to show that: . Substituting this result into equation (38a) it is found that P satisfies the constitutive relation: Formulas (42) and (44)  Hence, using again (41) we conclude that the magnetic permeability tensor must be symmetric: In summary, the constraints (41) reproduce well known results for reciprocal bianisotropic media, particularly they impose that the permittivity and permeability tensors are symmetric and that the magneto-electric coupling tensors satisfy ξ ζ = − EB EB T [42].

The number of plane wave modes and the ABC requirement
Let us consider a scattering problem wherein a plane wave propagating in a vacuum (or alternatively in a regular dielectric) illuminates a semi-infinite quadrupolar medium (region z > 0; the boundary is at the plane = z 0). It is well-known that because of translational invariance if the incoming wave varies in the x-and y-coordinates as e e k x k y i i being the propagation factor of a generic plane wave mode in the quadrupolar medium. Next, we determine the number of modes that can be excited in the quadrupolar medium, and link it with the number of boundary conditions at the interface.
The dispersion equation for the plane wave modes in a generic spatially dispersive material is [6,35]: If the nonlocal dielectric function ε ω k ( , ) ef is a polynomial function of k (as in equation (39), which will be assumed to hold in what follows) then P k ( ) D is a polynomial function of k.
Let us put where k k , x y 0 0 are determined by the excitation as previously discussed. We want to find the number of roots of ( ) . In case of spatial dispersion of second order, each element of the dyadic A in equation (47) is a polynomial of degree 2 or less in k z , and hence ( ) is necessarily a polynomial of degree 6 or less. Next, we prove that if the material response is such that Q zz is identically zero, i.e. if the bulk effective parameters are such that Q zz vanishes for any macroscopic electric field or electric field derivatives in the bulk region, then is a polynomial of degree 4 in k z . The condition that Q zz is identically zero for any field is equivalent to stating that ω ⋅ = q u k ( , ) 0 z z , which gives: To find the degree of ( ) we use equation (39) to obtain: where ε 0 is some dyadic independent of k z . From the reciprocity relations (41) it is possible to write: From here it can now be checked that when Q zz is identically zero in the bulk quadrupolar medium (equation (48)), the effective dielectric function is such that ε zz ef, is independent of k z and ε xz ef, , ε yz ef, , ε zx ef, , ε zy ef, are polynomials of degree 1 or less in k z . This implies that the dyadic A in equation (47) has a similar property and consequently, from the properties of the determinant operator, is a polynomial of degree 4. In summary, for quadrupolar media with weak spatial dispersion (equation (39)) the characteristic equation for the plane wave modes is generally a polynomial of degree 6.
However, when the material response is such that Q zz is identically zero in the bulk region (equation (48)), then ( ) becomes a polynomial of degree 4. Let us now discuss the implications of this result in the context of the scattering problem outlined in the beginning of this section. Physically, it is evident that for reciprocal media the number of physical channels associated with propagation along the +z-direction is the same as the number of channels associated with propagation along the −z-direction. Hence, the number of plane wave modes that propagate away from the interface is equal to the degree of ( ) P k z D0 divided by two. Thus, in general an incoming plane wave propagating in the dielectric can excite three distinct plane waves in the bulk quadrupolar medium. Because the reflected wave can be written in terms of two plane waves (there are two independent polarizations in a standard dielectric), we see that in general our scattering problem has five scalar unknowns (the complex amplitudes of the two reflected waves and the complex amplitudes of the three transmitted waves). To determine these five unknowns one needs five scalar boundary conditions. The vector boundary conditions (13)-which generalize the standard Maxwellian boundary conditions-only provide four scalar equations. Thus, the need for an additional scalar boundary condition (equation (32)) becomes manifest.
There is, however, an exception wherein an ABC is not required. If the material response is such that Q zz is identically zero in the bulk region (equation (48)) only two waves are excited in the quadrupolar medium, and hence the boundary conditions for the tangential fields (13) are sufficient, and the ABC should be redundant so that the boundary conditions are not over specified. Thus, if Q zz is identically zero in the bulk region, the parameter α Q in equation (32) must vanish to guarantee that the ABC is redundant.

Example of application
To illustrate the application of the developed theory, we consider a metamaterial formed by a periodic array of pairs of parallel metallic nanorods ( figure 2(a)). This metamaterial has been widely studied in the context of negative refraction at optical frequencies (e.g. [43][44][45]). Several works have highlighted that the nanorod pair metamaterial and other related designs have a strong electric quadrupolar response [15,24,45]. The origin of the quadrupolar response is well understood [24,45], and is related to the excitation of an anti-symmetric mode with currents in the two nanorods oriented in opposite directions ( figure 2(a)). From the geometry of the basic inclusion, it is expected that the anti-symmetric mode may be excited by a y-gradient of the xcomponent (parallel to the nanorods) of the electric field. In other words the anti-symmetric mode can be excited by a nonzero ∂ E y x . Moreover, from the definition (7b) it is seen that for the anti-symmetric mode the total quadrupole density can only have an yx-component. Thus, this discussion suggests thatˆ= ⊗ Q Q y x yx where Q yx responds to ∂ E y x . Below we report an analytical model that is consistent with this result. Note that because Q is neither antisymmetric nor symmetric, both an electric quadrupole density Q el and a magnetization vector must coexist ( figure 2(b)). The reason why ≠ Q 0 el is because the current loop associated with the nanowire pair is discontinued at the nanowire edges ( figure 2(b)).
In [45] an analytical model for the optical response of the bulk material was developed. Under the assumption that the metallic nanorods are oriented along the x-direction, are electrically thin, and are separated by a distance d along the y-direction, it was shown in [45] that the polarization vector, the electric quadrupole moment, and the magnetization vector in the metamaterial are such that≈ P P , , x z xy el, given by (it is assumed that ≪ k d /2 1 y ; our notations differ from those of [45]): x y x 0 where ω 0 , ω δ , γ, and B are parameters that depend on the dimensions, relative distance, material parameters and on the volume fraction of the nanorod pairs. Here, for simplicity, the values for these parameters are taken from the gold nanorod pair design reported in figure 2 and table I of [45], so that ω = a c / 2.78 0 , ω = δ a c / 1.55, γ = a c / 0.19, and = B 1.14, where = a 600 nm is period along the x-direction, and = d 65 nm. For more details the reader is referred to [45]. In figure 2(c) we show χ + and χ − , i.e. susceptibilities that determine the frequency dispersion of constitutive parameters of the metamaterial, as a function of the normalized frequency.
From equation (51a) it is seen that P x responds to both E x and ∂ E y x 2 . Thus, the polarization response is nonlocal. On the other hand, equation (8) and equation (51b) confirm that the strength of the electric quadrupolar response is comparable to that of the magnetization. Indeed, by combining the electric quadrupole with the magnetization as in equation (8), one obtains the total quadrupole moment density: The currents associated with the anti-symmetric resonance can be seen as a superposition of a magnetic dipole and an electric quadrupole. (c) Constitutive parameters for a particular metamaterial geometry as a function of frequency (following [45] Next, we characterize the plane waves in the bulk metamaterial. We restrict our attention to waves that propagate in the xoy plane with wave vectorˆ= + k k k x y x y , and with a magnetic field directed along the z-direction: . The most convenient way to obtain the dispersion of the photonic plane wave states is to use the spatially dispersive dielectric function ε ω k ( , ) ef , as discussed in section 4.4 (equation (47)). From  = − ⋅ Q P P g (equation (9)), equations (51a), (52) and ∂ = k i y y it is found that Therefore, the bulk material can be seen as a uniaxial spatially dispersive material characterized by the dielectric function ε ω k ( , ) ef . The reason why the only nonzero component of ε ef is the xx-component is because the wires are electrically thin, and hence interact weakly with a field oriented along the y-direction. In agreement with [45], it is found that the dispersion of the photonic states for wave propagation in the xoy plane with= The electric fieldˆ= ε ω Let us now consider the problem of wave scattering by a metamaterial slab with thickness L ( figure 3(a)). We assume that the incoming plane wave propagates in the xoy plane and is p-polarized= . We want to find the reflected and transmitted waves with the effective medium approach under the assumption that the macroscopic fields are finite in all space, so that the relevant boundary conditions are as in equation (18) (or equivalently, as in equations (13b), (19), and (20)). The results of effective medium theory are expected to be more accurate for small values of ωa c / and of k a inc . It is pertinent to admit that the optical axes of the metamaterial may be misaligned with the respect to the interface (the reason will be understood shortly). Thus, it is assumed that the metamaterial slab interfaces are normal to the ′ x -direction and parallel to the ′ y -and z-directions, such that the x-direction makes an angle φ with the ′ x -direction ( figure 3(a) The formulas in equation (57a) may be regarded as a generalization of the usual Maxwellian boundary conditions to quadrupolar media. On the other hand, equation (57b) is equivalent to = Q 0 nn and corresponds to the ABC (20). This ABC is not redundant unless φ = 0 or φ π = /2, i.e. unless the optical axes of the quadrupolar material are aligned with the interfaces. Thus, to illustrate the use of the ABC we need to suppose that the optical axes are misaligned with the interfaces, and in what follows the tilt angle is taken equal to φ =°45 . In such a case, equations (57) are equivalent to: represents the quadrupole moment evaluated at the metamaterial side of the pertinent interface.
The wave vector associated with the incident wave k inc can be written aŝ=    . In this manner, by imposing the continuity of ′ E y and H z at the two interfaces ′ = x 0 and ′ = x L and by imposing that Q yx vanishes at the metamaterial side of the two interfaces (equation (58)), it is possible to reduce the scattering problem to a 6 × 6 linear system with unknowns R T A , , i , (i = 1,…,4). We emphasize that the ABC = Q 0 nn is essential to solve this scattering problem, because otherwise the number of unknowns would be larger than the number of equations. Figure 3(b) shows the amplitudes of the calculated R and T for the case where the tilt angle is φ = 45°and = L a 2 2 . The solid lines correspond to incidence from left to right with θ i = 15°, and the dashed lines to incidence from left to right with θ = −°15 i . The dot-dashed blue line represents the transmission coefficient amplitude for θ = −°15 i when the quadrupolar resonance is ignored and χ − is set equal to zero. Note that in such a case an ABC is not required and the problem can be solved with standard methods. As expected, the effects of the quadrupolar response are negligible away from the quadrupolar resonance ωã c ( / 2). In this example the size of the unit cell at the quadrupolar resonance is of the order of λ a /3 0 and hence the effective medium theory is pushed to its limits of applicability.
It is also seen in figure 3 i . This is tantamount to saying that = S S 12 21 where S is the scattering parameter associated with ports 1 and 2 represented in figure 3(a). The symmetries of the S-parameters are a direct consequence of the Lorentz reciprocity theorem [46]. In the same manner, we numerically verified (not shown) that the transmission coefficient from right to left is exactly the same as the transmission coefficient from left to right = [ ] S S 13 31 . We also checked that when loss is removed (γ is set to zero) the boundary conditions ensure that + = R T 1 2 2 . These results provide a numerical check of our theory.

Conclusion
We investigated the boundary conditions for the macroscopic electromagnetic fields in quadrupolar media. Starting from first-principles microscopic considerations we derived the boundary conditions satisfied by the tangential components of the fields (equation (13)). Our results agree with other formulas previously reported in the literature derived using different methods [18,19]. We proved that the electromagnetic response of reciprocal quadrupolar media must satisfy a consistency condition (equation (28)). A solution for the consistency condition was found, and imposes that the bulk effective parameters have certain symmetries equation ((41)) and that the macroscopic fields satisfy an ABC at interfaces with regular dielectrics (equation (32)). The ABC depends on the bulk response of the involved materials and on a parameter (α Q ) which may depend on the particulars of the interface and internal degrees of freedom. We have demonstrated that if the macroscopic electromagnetic fields are sectionally continuous at the interface the parameter α Q must vanish and the ABC is simply = Q 0 nn (equation (20)). The need for an ABC was linked to the emergence of additional waves in the bulk quadrupolar media. To illustrate the application of the theory we studied the scattering of plane waves by a metamaterial slab formed by metallic nanorod pairs. Interestingly, when α ≠ 0 Q both the tangential electric field and the normal component of the Poynting vector may be discontinuous at the boundary (see appendix C). Even though such properties are rather exotic and it is tempting to adopt α = 0 Q as a universally valid result, it does not seem possible to completely exclude the possibility of having macroscopic fields with singularities, and only future numerical and experimental studies can clarify this important issue.