Polarizability Extraction for Waveguide-Fed Metasurfaces
Laura Pulido-Mancera Patrick T. Bowen Nathan Kundtz Mohammadreza F. Imani David Smith
physics.app-ph

APS/123-QED




 laura.pulido@duke.edu

Department of Electrical and Computer Engineering, Duke University, Durham, NC, 27708.

Kymeta Corporation, 12277 134th Court NE, Redmond, Washington 98052, USA

This work was supported by the Air Force Office of Scientific Research (AFOSR, Grant No. FA9550-12-1-0491)


 Department of Electrical and Computer Engineering, Duke University, Durham, NC, 27708.

 Kymeta Corporation, 12277 134th Court NE, Redmond, Washington 98052, USA



 Department of Electrical and Computer Engineering, Duke University, Durham, NC, 27708.



Kymeta Corporation, 12277 134th Court NE, Redmond, Washington 98052, USA



Department of Electrical and Computer Engineering, Duke University, Durham, NC, 27708.



We consider the design and modeling of metasurfaces that couple energy from guided waves to propagating wavefronts. This is a first step towards a comprehensive, multiscale modeling platform for metasurface antennas—large arrays of metamaterial elements embedded in a waveguide structure that radiates intro free-space—in which the detailed electromagnetic responses of metamaterial elements are replaced by polarizable dipoles. We present two methods to extract the effective polarizability of a metamaterial element embedded in a one- or two-dimensional waveguide. The first method invokes surface equivalence principles, averaging over the effective surface currents and charges within an element to obtain the effective dipole moments; the second method is based on computing the coefficients of the scattered waves within the waveguide, from which the effective polarizability can be inferred. We demonstrate these methods on several variants of waveguide-fed metasurface elements, finding excellent agreement between the two, as well as with analytical expressions derived for irises with simpler geometries. Extending the polarizability extraction technique to higher order multipoles, we confirm the validity of the dipole approximation for common metamaterial elements.  With the effective polarizabilities of the metamaterial elements accurately determined, the radiated fields generated by a metasurface antenna (inside and outside the antenna) can be found self-consistently by including the interactions between polarizable dipoles. The dipole description provides an alternative language and computational framework for engineering metasurface antennas, holograms, lenses, beam-forming arrays, and other electrically large, waveguide-fed metasurface structures.


Valid PACS appear here
                             

                              
Polarizability Extraction for Waveguide-Fed Metasurfaces
    David Smith
    July 14, 2017
========================================================







§  INTRODUCTION

The concept of metamaterials—artificially structured materials—has prompted the rapid development of new tools and techniques for controlling the propagation of waves. The metamaterial paradigm, in which an artificial medium is assembled from a collection of subwavelength elements with desired scattering characteristics, has had a profound impact across numerous scientific fields, including physics <cit.>, electromagnetic <cit.> and acoustic <cit.> wave phenomena, materials science, chemistry, engineering <cit.>, and nanoscience <cit.>. Metamaterials have provided a venue to tailor material properties in ways not feasible with conventional materials <cit.>, opening the door to unique and often exotic wave phenomena such as negative and near-zero refractive index materials  <cit.>. Furthermore, metamaterials research has lead to the demonstration of unprecedented devices, such as transformation optical structures and invisibility cloaks <cit.>, as well as superlenses<cit.>.  

The underlying metamaterial paradigm is that the behavior of waves propagating within a large (many wavelengths) composite medium can be understood from the properties of constituent elements—each subwavelength in dimensions—and their mutual interactions. The advantage of the metamaterial perspective is that the properties of each of the constituent elements can be determined exactly using a full-wave simulation over a relatively small, subwavelength, spatial domain. From these simulations, effective constitutive parameters can be retrieved, replacing the detailed current and field distributions within the small domain by just a few parameters such as the electric permittivity and magnetic permeability <cit.>. The wave propagation properties of the composite structure can then be modeled by solving Maxwell's equations directly, with effective constitutive parameters replacing the actual metamaterial structures. While the effective constitutive parameters obtained by numerical retrieval methods must be applied with considerable caution, retrieval methods have nevertheless been used with success in the design of many metamaterial structures  <cit.>. Replacing the details of an artificial medium with effective constitutive parameters facilitates device simulations and optimization cycles, vastly reducing the computational requirements since the individual elements are replaced by homogenized constitutive parameters. Complex metamaterial devices have been designed and demonstrated by this technique, including the transformation optical structures that rely on precise variations in material properties throughout a volume <cit.>.

Despite the compelling features of volumetric metamaterials and their unique and often unprecedented material properties, the applications for volumetric metamaterials have been limited. This limitation exists because some of the most intriguing properties of metamaterials occur near resonances in the individual elements—resonances that often impose bandwidth limitations and produce large resistive losses. Thus, waves propagating through any significant volume (even just a few wavelengths) of a metamaterial can be heavily attenuated. In addition, fabricating metamaterial elements to control electric and magnetic fields polarized in arbitrary directions, and assembling such elements throughout a volume, remains a major implementation challenge; typically, the properties of volumetric metamaterials have been demonstrated in highly constrained formats and proof-of-concept prototypes.

The difficulties associated with volumetric metamaterials are considerably reduced for metamaterial structures consisting of just a single layer or a few layers of metamaterial elements—also known as metasurfaces <cit.>. Being easier to design, model and implement <cit.>, metasurfaces have rapidly gained traction as a major subfield in metamaterials research <cit.>. As quasi-optical devices, metasurfaces provide control of reflection and transmission across the spectrum <cit.>, paving the way for advanced components such as flat lenses<cit.>, thin polarizers<cit.>, spatial or frequency filters <cit.>, and holographic and diffractive elements <cit.>. Used as coatings, metasurfaces can control the absorbance and emissivity of a surface, and thus have relevance to thermophotovoltaics <cit.>, detectors and sources <cit.>. Given the capabilities of metasurfaces to control waves, but without many of the limitations of volumetric metamaterials, metasurfaces have proven a good match for commercialization efforts, with many serious applications now being pursued, including satellite communications <cit.>, radar <cit.>, and microwave imaging <cit.>.

As with a volumetric metamaterial, the scattering properties of a metasurface can be characterized by a set of effective surface constitutive properties, which homogenize—or average over—the properties of many identical, discrete, metamaterial elements <cit.>. These effective medium properties relate to the discontinuity of the fields across the metasurface—approximated as having infinitesimal thickness—and are encapsulated in a set of generalized boundary conditions <cit.>. Effective surface constitutive properties are a convenient and natural description for free-standing metasurfaces illuminated by plane waves <cit.>; however, there are many contexts where a homogenization description is not the most convenient or the most accurate. For example, a metasurface interacts with a guided wave in ways not easily captured with a simple homogenized description. To better illustrate this point, consider a rectangular waveguide with one bounding surface patterned with a collection of metamaterial elements, each consisting of a pattern of voids or irises in the conducting surface. The volume of the waveguide just beneath the element behaves as if filled with a material possessing effective constitutive parameters introduced by the presence of the metamaterial element <cit.>. In addition, the waveguide-fed metamaterial elements also leak energy out of the waveguide, such that the entire structure can act as an aperture antenna <cit.>. Indeed, metasurfaces can provide nearly total control over the wavefront in a manner similar to phased arrays <cit.> and other aperture antennas, but often with advantages not available in other formats <cit.>. 

For metasurfaces, a weakness associated with homogenization techniques is that the element size and average spacing between metamaterial elements must be significantly subwavelength—a condition not necessarily satisfied in many situations. An element size of one-tenth to one-fifth of a wavelength is typical for many metamaterial structures, which implies the phase of the wave will have significant variation over the volume containing the element. Such metamaterials are said to exhibit spatial dispersion, and are properly characterized by constitutive parameters that depend on the wave vector in addition to the frequency, adding complication to the homogenization description. While numerical retrievals include spatial dispersion in the effective constitutive parameters, the retrieved parameters are specifically valid only in the exact arrangement simulated. For example, if the effective constitutive parameters are retrieved from a simulation of a cubic cell with periodic boundary conditions, the retrieved parameters will only be specifically valid for that medium, and not necessarily applicable when the same element is placed in a different context—for example, in a random arrangement.

To arrive at a more generally valid description of a metamaterial while still avoiding a full-wave simulation of the composite structure, we consider directly the properties of each metamaterial scattering element. The response of such a metamaterial element can generally be expressed in a series of induced electric and magnetic multipoles, typically dominated by the dipole term. The strength of the dipolar contribution is connected to an effective polarizability, which represent the coupling between the total dipole moment and the incident field on the element. To the extent that the higher order multipoles beyond the dipolar term can be neglected, the scattering from a collection of equivalent dipoles provides a near exact and computationally efficient model of the metamaterial structure. By assigning an effective polarizability to a metamaterial element rather than treating the metamaterial or metasurface as a continuous medium with constitutive parameters, it is possible to predict the overall response of the structure without any limitation on the element's periodicity or arrangement. The combination of polarizability extraction and the dipole representation forms an alternative, powerful modeling platform for metasurfaces and metamaterials. 


To obtain the polarizability of an arbitrary metamaterial element, a technique that can be applied is to first assume the element is part of an infinitely periodic medium <cit.>. The details of the periodic structure can then be replaced by periodic boundary conditions, so that the full-wave simulation domain extends only over a single cell of the structure. An effective polarizability of the element, which includes the contributions from all other elements in the infinite array, can then be extracted from the computed field or charge/current distributions. Finally, the intrinsic polarizability can be determined using the Lorentz formula that relates the intrinsic and effective polarizabilities <cit.>. 

Our aim here is to extend the dipole model as an analytical tool for waveguide-fed metasurfaces. Unlike the free-standing metasurface or volumetric metamaterial, for which each metamaterial element can be reduced to a free space dipole, an individual metamaterial element patterned in a waveguide also interacts with the waveguide structure. Even if the metamaterial element is reduced to a single free space dipole, the image dipoles induced within the waveguide structure must then be taken into account to obtain an accurate expression for the waveguide-fed metamaterial element <cit.>. Rather than working through this complication, we can instead apply a numerical polarization extraction procedure using the waveguide modes, arriving at an effective polarizability that includes all of the waveguide interactions. Once having established this dipole description of a waveguide-based metamaterial element, we can subsequently find the properties of the composite waveguide-fed metasurface accurately, efficiently and quickly, for example using the dipole model previously described.

We approach this problem as follows: In Section <ref> we introduce the polarizability framework for a metamaterial element embedded in a waveguide and summarize the self-consistent dipole model for metasurfaces. In Section <ref> we apply the equivalence principle to the waveguide-fed metamaterial, and derive integrals relating the equivalent current densities to the effective dipole moments. The results of the direct integration method are used as the basis for comparison with the second method we present in Section <ref>, in which the polarizability is obtained using the scattering (S) parameters of the element when placed in a waveguide. Both methods have been used in the context of numerical retrieval of effective constitutive parameters for volumetric metamaterials, with the former method related to field averaging <cit.>, while the latter method related to the well-known S-parameters retrieval method <cit.>. In Section <ref> we perform polarizability extractions for different waveguide-fed metasurface geometries: circular iris <cit.>, elliptical iris <cit.>, iris-coupled patch antenna <cit.> and the complementary electric inductive-capacitive metamaterial resonator(cELC) <cit.>. We extend the polarizability extraction method to two-dimensional (2D-) waveguide structures in section (<ref>), where the S-parameters cannot be used to retrieve polarizability. To address this problem, we outline an alternative method based on the mode expansion of cylindrical waves propagating through the waveguide. This framework is particularly advantageous for modeling and designing planar structures <cit.>. We also demonstrate that the extracted polarizability of a metamaterial element depends on the geometry of the waveguide in which it is embedded. As a means of confirming the accuracy of the dipole approximation, we express the induced fields from a metamaterial element as a multipole expansion, comparing the relative strengths of the expansion coefficients. The result of this analysis shows that, indeed, the dipole term dominates the response, justifying the dipolar description and use of the dipole model. We conclude by examining the potential application of the polarizability extraction for different metasurface-based devices.



§ SELF-CONSISTENT DIPOLE MODEL FOR METASURFACES 


Our conceptual picture of a waveguide-fed metasurface is that of a collection of polarizable dipoles, each of which accounts for the scattering associated with an actual metamaterial element. Using the polarization extraction technique to be presented, a metamaterial element is reduced to an effective electric dipole moment and an effective magnetic dipole moment, 𝐩 and 𝐦, which are proportional to the local electric or magnetic field at the center of the element multiplied by a coupling coefficient, termed the dynamic polarizability <cit.>. It is referred to as “dynamic" since it describes the element's response due to a time-varying incident field. A time-dependent electric field can induce solenoidal currents and thus can give rise to a magnetic polarization in addition to the electric polarization. To properly take into account the co- and cross-coupling of the electromagnetic fields excited within the metamaterial element, the dynamic polarizability should be represented as a tensor, so that 𝐩= α̅̅̅_ee𝐄^loc and  𝐦 = α̅̅̅_mm𝐇^loc, where 𝐇^loc and 𝐄^loc are the local electric and magnetic fields at the center of the element.  

It is worth noting that the polarizability of a metamaterial element is not an inherent property, but rather depends on the environment, and therefore may be considered a nonlocal property. To better illustrate this point, consider a single passive magnetic dipole m placed at position r_0 in an environment that is described by the Green's function G(r,r'). This configuration is illuminated by an incident magnetic field H^0. The total field H^tot everywhere throughout the environment is thus

    H^tot=H^0+G(r,r_0)m.


There are two possible ways to define the polarizability of a dipole <cit.>. One is the effective polarizability, defined in terms of the incident wave as m=α_mH^0, and the other is the inherent polarizability, defined in terms of the total field m=α̃_mH^tot. The latter case would be a difficult definition to uphold consistently for a point particle, because the real part of the Green's function G(r,r_0) diverges as r→r_0. (Some authors have developed ways of dealing with the singularity and relate it to the radiation of the dipole <cit.>.) 

The imaginary part of the Green's function at the origin, however, is well-known to have an important physical meaning related to the radiation reaction force on the dipole <cit.>. Defining the polarizability as m=α_mH^0, we consider the average power absorbed by the dipole according to Poynting's theorem: 

    P_abs=(1/2)∫Re{J_m^*·H}dV

where J_m=iωμ_0mδ^(3)(r-r_0) is the magnetic current associated with the magnetic dipole,  ω is the angular frequency, ^* represents complex conjugate, and μ_0 is the permeability of free space. If there are no Ohmic losses in the system, then all power scattered by the dipole is re-radiated, and consequently P_abs=0. Since the magnetic field H that appears in Eq. <ref> is the total magnetic field given in Eq. <ref>, P_abs=0 implies that

    Im{m^*·H^0-m^*G(r,r_0)m}=0.

Substituting m=α_mH^0 and assuming that the polarizability is isotropic, we obtain a relationship between the imaginary part and the real part of the polarizability as

    Im{α_m}=|α_m|^2Im{m^*G(r_0,r_0)m/|m|^2},

which must be satisfied in order for the conservation of energy to be upheld. Equation (<ref>) implies that Im{1/α_m}=Im{m^*G(r_0,r_0)m/|m|^2}. Therefore, after some simple algebraic manipulation, Eq. <ref> can be recast as  

    α_m=α̃_m/1+iα̃_mIm{m^*G(r_0,r_0)m/|m|^2}.

Here, α̃_m is the inherent polarizability of the dipole, whereas α_m, which is defined as the ratio between the dipole moment and the incident field, is dependent on the local environment. If a dipole is placed in free space, then a Taylor series expansion of the Green's function shows that Im{G(r_0,r_0)}=k^3/6π, and this yields the radiation reaction to the polarizability of a dipole in free space

    α_m=α̃_m/1+iα̃_mk^3/6π.


The expression in Eq. <ref> is often known in the literature as the radiation reaction correction or the Sipe-Krankendonk relation <cit.>. If instead the dipole is placed just above an infinite ground plane, the dipole radiates twice as much energy, and so Im{m̂G(r_0,r_0)m̂}=k^3/3π. If the dipole is a complementary metamaterial element embedded in a waveguide wall, then it will radiate both into the upper half space and into the waveguide, and so the radiation reaction correction would need to take into account both scattered fields. Considering that this correction must account for half of the radiation in free-space, as in Eq. <ref>, and half of the radiation inside the waveguide, as derived in appendix <ref>, the corrected polarizability has the form


    α_m=α̃_m/1+iα̃_m (k^3/3π + k/ab).


Equation (<ref>) has significant implications on any polarizability extraction method that deals with waveguide integrated metamaterial elements. For example, if the static polarizability of an element is calculated using Bethe theory (as discussed in section <ref>), then the radiation reaction correction will be different depending on whether that element is placed in a 2D waveguide or a cavity, and so the proper correction as in Eq. <ref> will need to be applied in each environment. 

Once the polarizabilities are determined, an accurate, self-consistent description of the scattering from the collection of dipoles can be obtained using the dipole model–including the Green's function associated with the specific environment–which expresses the relationship between the local fields 𝐄^loc(𝐫_i) in terms of the incident fields 𝐄^0(𝐫_i) and the fields scattered from all of the dipoles. These equations can be written as

    𝐄^loc(𝐫_i) =𝐄^0(𝐫_i) + 
    ∑_j≠ i𝐆_ee(𝐫_i,𝐫_j) 𝐩(𝐫_j) + 𝐆_em(𝐫_i,𝐫_j) 𝐦(𝐫_j)

where 𝐆_ee, 𝐆_em correspond to the electric components of the dyadic Green's function. An equivalent expression can be derived for the magnetic field. Examining the expression in Eq. <ref>, it can be seen that these coupled equations capture the interaction of the incident wave with each of the meta-atom (through the j=0 terms) as well as interaction between different elements (the summation term). Note that we do not solve these equations in the present work, but include them to provide the complete modeling framework. The development of the waveguide Green's function and simulation of various waveguide-fed metasurfaces will be presented separately.



§ POLARIZABILITY EXTRACTION IN A RECTANGULAR WAVEGUIDE: DIRECT INTEGRATION


We start by considering an arbitrarily shaped iris etched into the upper conducting surface of a rectangular waveguide, as shown in Fig.<ref>. The coordinate system is chosen so that the propagation direction is in the z-direction; a corresponds to the width of the waveguide along the x-axis, and b corresponds to the height along the y-axis. The guided wave couples to the iris, which radiates a portion of the incident wave into the free space region. A common methodology to solve the radiated field by such a structure is the surface equivalence principle. This principle states that the electric field on the boundary of a domain can be represented as a magnetic surface current K_m=E×n̂, while the magnetic field on the boundary can be represented as an electric surface current K_e=n̂×H, where 𝐄 and 𝐇 are the total fields on the surface of the domain, and 𝐧̂ is the normal to the surface. Using K_m and K_e and the corresponding Green's functions, one can determine the field within the domain. 

Applying this principle to the geometry of Fig. <ref>, we see that the tangential electric field is zero everywhere on the waveguide surface except over the void regions defining the iris; if the iris is deeply subwavelength, then the field scattered into the far-field may be approximated by just the first term of the multipole expansions of K_e and K_m. Hence, the dipole moments representing the iris can be calculated as <cit.>


    𝐩 = ϵ_0 n̂∫r·𝐄^tan da 
    𝐦 = 1/i μω∫n̂×𝐄^tan da .


The integration is performed over the surface of the iris, n̂=ŷ is the vector normal to the top surface, and 𝐄_tan corresponds to the tangential field at the surface of the iris. It is worth noting the tangential magnetic field is not zero over the surface of waveguide; however, the tangential magnetic field corresponds to an electric current density parallel to a metallic wall, and by image theory its effect can be ignored. 


Since the effective dipole moments are proportional to the incident fields 𝐄^0,𝐇^0, we can define the effective polarizabilities as







    𝐩 = ϵ_0  α̅̅̅^e_eff𝐄^0  𝐦 = α̅̅̅^m_eff𝐇^0.


In the most general case, each iris can be described by an electric polarizability tensor α̅̅̅^e_eff and a magnetic polarizability tensor α̅̅̅^m_eff. These tensors are symmetric and would thus normally have six unknown parameters, each of which would need to be extracted from a full-wave simulation. In free space, these components can be found by computing the scattered fields in all directions. In the waveguide geometry considered here, such simulations are not possible. Instead, we take advantage of the geometry of the iris and the symmetry of the excitation field into account. For example, when the symmetry axes of iris are coincident with the symmetry axes of the waveguide, we can assume that the polarizability tensor is diagonal. Further, the boundary conditions of the waveguide require that the incident tangential electric field E_x^0 and E_z^0 and the normal component of the magnetic field H_y^0 be equal to zero on the surface of the waveguide. Therefore, the tensor components α_px, α_pz, and α_my can never be excited, and hence we can assume them to be zero. The polarizability tensors then reduce to α̅̅̅^e_eff=diag(0,α_py,0) and α̅̅̅^m_eff=diag(α_mx,0,α_mz). Hence, the polarizability extraction is simplified to finding three unknowns (α_py, α_mx, α_mz). One further unknown can be removed if the iris is placed where the magnetic field of the incident mode has a null in the z-component. In that case, only α_py and α_mx are relevant to the problem. By using Eq. <ref> and Eq. <ref> and the previously described assumptions, the characteristic polarizabilities are


    α_ey = ∫∫(xE_x + zE_z)dxdz/E^0_y
    α_mx = 1/i μω H^0_x∫∫ E_z dxdz .


Equations Eq. <ref> provide a simple method to calculate the polarizability of an element from a full-wave simulation of the fields of the element embedded in a waveguide structure.



§ POLARIZABILITY EXTRACTION IN A RECTANGULAR WAVEGUIDE: SCATTERING PARAMETERS


While calculating the polarizability of a metamaterial element by means of Eq. <ref> provides a physically accurate characterization, the integration over the surface of the element can be cumbersome to perform for all desired frequency points and for arbitrary geometries. In many instances, this integration may also be subject to numerical inaccuracies due to singularities near edges or coarse meshing, as it is especially the case for resonant elements such as those examined in section <ref>. Instead of the direct integration, in this section we consider the extraction of the polarizabilities from the fields scattered by the element into the waveguide. For this calculation, we apply coupled mode theory to determine the coupling of the element embedded in a waveguide to the forward and backward scattered fields within the waveguide. The fields inside the waveguide at any plane of constant z (along the propagation direction) can be expanded as a discrete sum of orthogonal modes. These modes are defined as  <cit.>


    E_n^+   =(E_n t(x,y)+E_n z(x,y))e^-iβ_n z
    H_n^+   =(H_n t(x,y)+H_n z(x,y))e^-iβ_n z
    E_n^-   =(E_n t(x,y)-E_n z(x,y))e^iβ_n z
    H_n^-   =(-H_n t(x,y)+H_n z(x,y))e^iβ_n z


where E^-_n and E^+_n are respectively the waveguide modes traveling in the backwards and forwards directions. The subscript “t" refers to the component of the fields that are transverse to the direction of propagation, and β_n is the propagation constant of the nth mode. The mode normalization used in Eq. <ref> is defined from the integral over the cross section of the waveguide, such that 

    ∫E_n·E_mda=δ_mn,

where δ_mn is 1 for n=m and 0 otherwise. Note that the electric fields of the modes are dimensionless, and the magnetic fields have units of inverse impedance as

    ∫H_n·H_mda=δ_mn/Z_n^2,

where the wave impedance Z_n is defined as a normalization constant for each  mode as

    Z_n=1/∫E_n×H_n·n̂ da.


In Eq. <ref> the integration is over the cross sectional surface of the waveguide, i.e. the surface representing Port 1 in Fig. <ref>. 

Consider a metamaterial element placed at the center of the top plate of the waveguide. We assume that the incident field is the forward-propagating fundamental mode—coming from Port 1—with unit amplitude E^0+, as shown in Fig. <ref>. When the metamaterial element is present, it couples and scatters to all modes. While the element has a finite size, for points inside the waveguide that are few wavelengths away, the element is well-approximated as a point scatterer placed at z=0 (location of the element). In the absence of the metamaterial element, the total field is simply the incident field, which is identical in both the forward and backward directions. As a result, we express the modal decomposition of the total (both incident and scattered) fields into backwards propagating modes at Port 1 (z=-Δ) as

    E^-=E_0^++∑_n A_n^-E_n^-,

where A^-_n are the amplitudes of the modes scattered by the element in the backwards direction, and n is the mode number. Similarly, a modal decomposition of the fields in the plane of z=+Δ into forward propagating modes yields

    E|^+   =E_0^++∑_n A_n^+E_n^+

where A^+_n are likewise the mode amplitude coefficients of the scattered field by the element in the forward direction, and the incident field E_0^+ has been written as a separate term. In these calculations, Δ can be any distance as long as it is larger than the size of metamaterial element.



We first consider a volume within the waveguide that encompasses the metamaterial element (bounded by the two ports). The incident field impinging on the element will induce a set of fields that we denote as E and H. Within the coupled mode formulation, these fields can be related to the waveguide modes through Poynting's theorem, or 

    ∇· (𝐄×𝐇^(±)_n - 𝐄^(±)_n ×𝐇)= J_e·E^(±)_n-J_m·H^(±)_n.

We integrate Eq. <ref> over the volume V the portion of the waveguide between the two ports, and applying the divergence theorem, Eq. <ref> becomes

    ∫_S (𝐄×𝐇^(±)_n - 𝐄^(±)_n ×𝐇)·𝐧̂ da = 
       ∫_V J_e·E^±_n-J_m·H^±_n dV

where S is the closed surface that encloses V and 𝐧̂ is an outwardly directed normal. Since the waveguide walls are assumed to be perfectly conducting, the only nonzero contributions to the surface integrals arise from the surfaces representing Port 1 and Port 2 (depicted in Fig.<ref>), and the surface of the metamaterial element. Since the field E can be written as a fictitious magnetic surface current through K_m=E×n̂ and H can be related to a fictitious electric surface current in the same way, then the surface integral in Eq. <ref> indicates the manner in which the effective dipoles representing the metamaterial element couple to each of the waveguide modes, as might be expected from Lorentz reciprocity. 

To obtain the amplitude coefficients A_n^(±), we assume there are no current sources in the volume, implying that the volume integral in Eq. <ref> vanishes. Substituting the expansions of the fields in Eq. <ref> and Eq. <ref> into Eq. <ref> and using the orthogonality relations in Eq. <ref>, we obtain the scattered amplitude coefficients as an overlap integral of the waveguide mode fields with the total field taken over the surface of the metamaterial element. More explicitly, the amplitude coefficients can be found as 

    A_n^(±)= Z_n/2∫_element (𝐄×𝐇^(∓)_n-𝐄^∓_n ×𝐇)·𝐧 da.
 
In an alternative approach, the electric field in the aperture could be considered zero and replaced by an equivalent electric and magnetic surface current, according to the equivalence principle. In this case, the surface integral vanishes everywhere except over the surfaces of the ports, but the volume integral over the meta-atom becomes a surface integral of the equivalent surface currents K_m=E×n̂ and K_e=-H×n̂. Using Eq. <ref> and Eq. <ref> and invoking orthogonality, we obtain

    A_n^(±)= Z_n/2∫_element (𝐊_e ·𝐄^(∓)_n-𝐊_m ·𝐇^(∓)_n)da.
 
Since the metamaterial element is deeply subwavelength, the fields of the waveguide modes can be expanded in a Taylor series around the center of the element. The lowest order term is constant over the surface of the element, yielding

    A_n^(±)= Z_n/2[   𝐄^(∓)_n(r_0)·∫_element𝐊_e da
       -𝐇^(∓)_n(r_0)·∫_element𝐊_m da].
 
As previously stated in Eq. <ref>, the two integrals in Eq. <ref> are proportional to the electric and magnetic dipole moments, p and m. Therefore, the final expression for the amplitude coefficients in terms of these dipole moments is given by


    A_n^+=iω Z_n/2(E^-_n·p - μ_0 H^-_n ·m) 
    
    A_n^-=i ω Z_n/2(E^+_n·p - μ_0 H^+_n·m).


Equation <ref> shows that the dipole moments are related to the incident fields, which in turn can be expanded in terms of eigenmodes. Since the incident field is the fundamental mode, the polarizability is defined by


    p=ϵα̅̅̅^eff_eE_0^+ 
    m=α̅̅̅^eff_mH_0^+.


Due to the symmetry of the fields in the rectangular waveguide, the α_mz component cannot be excited by the z-component of the magnetic field. Hence, Eq. <ref> reduces to two coupled equations with two unknowns: α_py and α_mx, which can be recast as


    A_n^+=iω Z_n/2(ϵα_pyE^+_0yE^-_ny - μ_0 α_mxH^+_0xH^-_nx) 
    
    A_n^-=i ω Z_n/2(ϵα_pyE^+_0yE^+_ny - μ_0 α_mx H^+_0xH^+_nx)


Considering the orthogonality of the eigenmodes and the symmetry properties of the electromagnetic fields—the transverse components of the electric field are symmetric under a flip of direction (i.e. E_ny^-=E^+_ny), while the magnetic field is antisymmetric (i.e. H_nx^-=-H^+_nx)—we can solve Eq. <ref> in order to find the polarizabilities as


    α_py= 2/iω Z_n(A_0^+ + A_0^-)/ϵ (E^+_0y)^2
    α_mx= 2/iω Z_n(A_0^+ - A_0^-)/μ (H^+_0x)^2.


For the fundamental mode, the normalized fields and impedance at the dipole location are given by

    |E^+_0y|^2= 4/ab   |H^+_0x|^2= 4β_10^2/ab Z_0^2 k^2   Z_0 = η k/β_10

where η is the vacuum impedance. Furthermore, the amplitude coefficients A^+_0 and A^-_0 correspond to the amplitude terms for the fundamental mode of the scattered fields in the forward and backward directions. Therefore they are directly related to the scattering parameters with respect to each port: the reflected field, related to A^-_0 is proportional to the reflection coefficient i.e. S_11, while the transmitted field related to A^+_0 is proportional to the transmission coefficient S_21 and the incident field in the forward direction. More explicitly, these relationships are expressed by

    A^-_0= S_11   A^+_0= S_21-1.

Taking into account Eq. <ref> in conjunction with Eq. <ref> and Eq. <ref> it is possible to find the final expression for the polarizabilities as 


    α_py   =-iabβ_10/2k^2(A_0^++A_0^-) = -iabβ_10/2k^2(S_21+S_11-1) 
    α_mx   =-iab/2β_10(A_0^+-A_0^-) = -iab/2β_10(S_21-S_11-1) .


Equation <ref> provides the polarizabilities of any metamaterial element embedded in a rectangular waveguide in terms of the scattering parameters, which can be obtained from direct measurement or full-wave simulation. However, it is important to note that the scattering parameters must be de-embeded to the plane of the metamaterial element before using them in Eq. <ref>. This is a straightforward process well-known in literature <cit.>. Another important point to note is that we have assumed ports which only excite/represent single mode. This condition should be applied when simulating these structures in numerical solvers. More importantly, since it is cumbersome in experiment to excite purely the fundamental mode, the ports should be placed at least a few wavelengths in distance away from the metamaterial element to ensure the non-propagating higher order modes have decayed. 

The equations in Eq. <ref> are similar to the expressions found for the effective polarizabilities of metamaterial elements in periodic metasurfaces <cit.>. This relationship means that a single element in a rectangular waveguide acts as a dipole whose response is equivalent to the response of the element in a periodic metasurface. The details of this equivalence will be discussed in a future paper.



§ SIMULATED RESULTS FOR THE RECTANGULAR WAVEGUIDE


Using a full-wave simulation, we can extract the polarizability of arbitrary metamaterial elements patterned into rectangular waveguides from the two different approaches described in sections <ref> and <ref>. For both extraction techniques, a single full-wave simulation in CST Microwave Studio is performed assuming a waveguide designed to operate over frequencies in the X-band (8-12 GHz). The waveguide dimensions are a=21.94  mm, b=5  mm, L=22.7  mm and thickness 1.27  mm. We perform this simulation for several different metamaterial element geometries: circular iris, elliptical iris, iris-coupled patch antenna, and the cELC resonator. 

In addition to the methods described above, the dipole moments of simple geometries, such as an elliptical iris may be also obtained from the static dipole moments of general ellipsoidal dielectric and permeable magnetic bodies. Consider an elliptically shaped aperture with the major axis along the x-direction and minor axis along the z-direction. Let the major radius be l_1 and minor radius l_2. In the static limit, the polarizabilities of such an elliptical iris are given by <cit.>
 
 
    α̃_̃m̃x̃ = 4π l_1^3 e^2/3[E(e) - K(e)]
    α̃_̃m̃z̃ = 4π l_1^3 e^2(1-e^2)/3[E(e) - (1-e^2)K(e)]
    α̃_̃ẽỹ= -4π l_1^3 (1-e^2)/3E(e),

 
where e= √(1-(l_2/l_1)^2) (assuming l_1>l_2) is the eccentricity of the ellipse, and K(e) and E(e) are the complete elliptic integrals of the first and second kind, respectively. If e=0 these expressions reduce to the static polarizabilities of circular irises. Using these static expressions for the dynamic polarizability will violate conservation of energy since it lacks the radiation damping term of dynamic polarizability. These expressions may be corrected by the radiation term, as described in Eq. <ref>, assuming the real part given in Eq. <ref>. 

Figure <ref> shows the polarizability of simple circular and elliptical irises computed using the two methods described in this paper as well as the theoretical methods of Eq. <ref> and Eq. <ref>. As shown, excellent agreement between the analytical expressions and the numerical extractions is obtained, verifying the proposed methods. Since the circular iris considered here does not possess a resonance, it is expected that the dynamic polarizabilities extracted numerically  are well-approximated by the theoretical expressions. Next, we examine the case of an elliptical iris, as shown in Fig.<ref>b. The extracted polarizabilities computed using the two numerical extraction methods of previous section exhibit excellent agreement. However, the analytical exhibits significant deviation from the numerical methods. This is expected since the elliptical iris supports a resonance over the frequency band of interest, which is not captured in the analytical expressions derived for the static field. This case further highlights the need for a precise numerical method to compute the polarizability of a metamaterial element.


It is worth noting that an elliptical iris etched in a rectangular waveguide is well-known in the antenna engineering community as the unit cell of a slotted waveguide antenna (SWA). SWAs are particularly attractive due to their advantages in terms of design simplicity, weight, volume, power handling, directivity, and efficiency <cit.>. SWAs rely on the gradual leakage of the guided mode through the slots. The metasurfaces considered throughout this paper also share this feature, since they also leak energy from the guided wave through the metamaterial elements. In other words, the methodology developed in this paper can also be applied to leaky wave antennas for modeling, simulation, and designing <cit.>.

While the circular and the elliptical irises may be analyzed using analytical expressions for the polarizabilities derived in the static limit, such closed-form expressions are not available for most metamaterial designs. For example, an element of potential interest in the design of metasurface antennas is the iris-fed patch, shown in Fig.<ref>a <cit.>. The inclusion of the metallic patch above the iris enhances the resonant response of the element, as exemplified by the narrower and stronger resonant response. Another common metamaterial element is the cELC, shown in Fig.<ref>b, commonly used in metasurface antenna designs. The resonant response of the cELC is highly susceptible to variations in its geometry <cit.>. For both elements, we observe excellent agreement between the two numerical methods of equation Eq. <ref> and equation Eq. <ref>, as shown in Fig. <ref> and <ref>.



The geometry of metamaterial elements can be quite complicated, such that the numerical integrals required in Eq. <ref> are likely to yield inaccuracies. For this reason, the extraction based on computing the waveguide scattering parameters is likely to be more reliable and easier to implement; moreover, this method can also be used in measurements on fabricated samples. It worth noting that in all of the results presented in Fig.<ref> and <ref>, it can be seen the electric polarizability is much smaller than the magnetic polarizability—in fact, three orders of magnitude smaller. This phenomenon is expected considering that the geometry under study corresponds to a small opening in a metallic wall. 





§  POLARIZABILITY EXTRACTION IN A PARALLEL PLATE WAVEGUIDE


In this section, we consider the case of a 2D waveguide and examine a metamaterial element that is etched on one wall of a parallel plate waveguide. While the polarizability extraction method based on the equivalence principle holds for the 2D waveguide-fed element, the nature of the waveguide modes changes substantially from the formulation in section <ref>. In this section, we modify this extraction technique to be applicable to planar waveguide systems. As previously described in section <ref>, a metamaterial element scattering into a waveguide can be described in terms of a sum of waveguide modes. Because the element is placed in the upper surface of the waveguide, the boundary condition dictates the tangential electric field and the normal magnetic field to be zero and the element can only couple to the transverse magnetic (TM) modes. Since the natural symmetry of the system is cylindrical, mode decomposition is simpler if we use cylindrical coordinates (r,θ). 

Setting the origin of the coordinate system to the center of the metamaterial element, the z-components of the scattered electric field—for the TM modes characterized by the (m,n) indices— are given by


    E^sc_z,c=β_m/k H_n^(2)(β_m r)cos(nθ) 
       E^sc_z,s=β_m/k H_n^(2)(β_m r)sin(nθ)


where the subscripts “c" and “s" refer to modes that have angular dependence cos(nθ) and sin(nθ), respectively. The propagation constant is given by β_m=√(k^2-(mπ/h)^2), where h is the height of the waveguide. Invoking the superposition principle, the total solution for the z-component of the scattered electric field can be expressed as


    E_z=∑_n∑_mA_mn^s E_z,s^mn+A_mn^c E_z,c^mn

When h<π/k, only the m=0 mode is propagating, and in this case the electric field at all points where r≫ h/π is dominated by the m=0 mode. Therefore we can reduce Eq. <ref> to 

    E_z=∑_nA_n^s E_z,s^0n+A_n^c E_z,c^0n.

The m=0 modes are given by 
 

    E_z,c^0n=H_n^(2)(kr)cos(nθ) 
       E_z,s^0n=H_n^(2)(kr)sin(nθ).



The amplitude coefficients, A_n, can be found from the scattered electric field E_z using the orthogonality of the {sin(θ),cos(θ)} basis. By integrating over a circle of radius r centered at the origin of the metamaterial element, as shown in Fig.<ref>, it is possible to define the amplitude coefficients as



    A_n^s=lim_r →∞1/π H_n^2(kr)∫_0^2π E_z(r,θ)sin(nθ) dθ
       A_n^c=lim_r →∞1/π(1+δ_n0)H_n^2(kr)∫_0^2π E_z(r,θ)cos(nθ) dθ.




To better illustrate the utility of equations Eq. <ref>, we consider a lossless parallel plate waveguide, fed by a cylindrical source oriented in the z-direction,, as shown in Fig. <ref>. The source is placed far enough from the metamaterial element to avoid evanescent coupling. A single cELC resonator is assumed to be patterned on the center of the top plate of the waveguide, with the same geometrical parameters as the cELC presented in the previous section. Since the full-wave simulation domain represents the total field instead of the scattered field,  this structure is simulated with and without the cELC, and the difference of the two simulation results are taken to obtain the scattered field due to the metamaterial element, such that E_z^sc= E_z^tot -E_z^0 at the plane z=h/2. Once the scattered field is computed, the integration outlined in Eq. <ref> is performed to find the amplitude coefficients. The integration radius is selected electrically large enough so that the evanescent modes have decayed —it is also ensured the integration curve does not contain the cylindrical source. 


Figure <ref> shows the magnitude of the amplitude coefficients for the scattered fields computed for the meta-atom shown in Fig.<ref>. To better illustrate the physics behind these coefficients, we apply Poynting's theorem Eq. <ref> which directly links the amplitude coefficients and the dominant dipole moments of the metamaterial element. In contrast to the 1D case examined in section <ref>, the amplitude coefficients are not related to the scattering parameters, but rather to the scattered fields, by means of Eq. <ref>. Moreover, while the location of the metamaterial element in the rectangular waveguide limits the calculation of a single component of the magnetic polarizability, such limitation disappears in the case of the planar waveguide. 

The cylindrical wave propagating through the waveguide may excite the two tangential components of the magnetic polarizability, which leads to a more complete characterization of the polarizability tensor of the element. The scattered fields generated by the metamaterial element can be represented as the sum of the moments of the surface current J^n_m multiplied by the different eigenmodes of the scattered fields shown in Eq. <ref>, more explicitly, this relationship is given by 

    E_z=m_xZ_0k^2/4h E_z,s^01+m_yZ_0k^2/4h E_z,c^01+-ip_zk^2/4hϵ_0 E_z,c^00

A direct mapping between Eq.<ref> and Eq. <ref> demonstrates that the first three amplitude coefficients, {A_0^c, A_1^c, A_2^s} are directly related to the three dominant dipole moments as <cit.>


    m_x=A_1^s4h/Z_0 k^2   m_y=A_1^c4h/Z_0 k^2   p_z=A_0^ci4hϵ_0/k^2.


As shown in Fig. <ref>, the predominant amplitude mode is A^s_1, which is directly associated with m_x, while the amplitude of the modes A^c_0 and A^c_1, associated with p_z and m_y, are significantly smaller—by two orders of magnitude. The effective polarizabilities given the incident wave due to the line source, can be directly obtained from their corresponding dipole moments given in Eq. <ref> as








    α_z^p=p_z/E_z^0  α_xy^m=m_y/H_x^0  α_xx^m=m_x/H_x^0.


For clarification, the double indices on the polarizabilities represent the entry in the polarizability tensor; for example, α_xy represents the component of the polarizability that generates a dipole moment oriented in y, due to the x component of the incident field. In order to find all three components of this tensor, it is necessary to rotate the cELC by π/2, and perform the same extraction technique. The electric polarizability and two of the components of the magnetic polarizability tensor are thereby obtained and shown in Fig.<ref>a. As shown, excellent agreement for the magnetic polarizabilities is obtained between the two numerical polarizability extraction methods. For this particular example, note that the polarizability α^m_yy has a resonant response out of the X-band, but its geometry can be modified such that it has both resonances in the same band <cit.>. In the case of the electric polarizability (Fig.<ref>b), the numerical values obtained are significantly smaller, which makes it susceptible to numerical inaccuracies when the integration in Eq. <ref> is performed. 



It is important to highlight in this example that the term A^s_3 is associated with the quadrupole moment. This term has been traditionally neglected in most metamaterial design strategy. The example at hand provides a useful framework to examine the contribution of the quadrupole term and the error introduced by neglecting it. To study the impact of the quadrupole term, we compared the simulated scattered field (shown in Fig. <ref> first row) with its theoretical expression Eq. <ref> up to only the dipolar contribution, as shown in the second row of Fig.<ref>. We observe excellent agreement between the two rows, confirming the assumption that the main contribution of the scattered field is dipolar. The physical implications of this result can be understood by calculating the difference between the scattered fields from a full-wave simulation and from the analytical expression in (Eq. <ref>). As shown in the third row of Fig. <ref>, the error due to assuming the dominant dipolar term is several orders of magnitude smaller than the amplitude of the scattered field, and the largest discrepancy is observed within the close vicinity of the metamaterial element. This result, in conjunction with the amplitude coefficients shown in Fig.<ref> also demonstrates that most of the radiation is associated with the dipolar term and higher order modes can be ignored. However, if the elements are placed at distances where these higher order modes have not decayed, they can alter the coupling between the two meta atoms and change the total scattered fields inside the waveguide. 



Another interesting point to highlight is that while the geometric characteristics of the cELC used in Fig.<ref>c and <ref> are the same, its resonant response as manifested by the effective polarizability changes depending on the host waveguide, as can be observed by comparing Fig.<ref>d and Fig.<ref>a. Therefore, the electromagnetic response of a metamaterial element embedded in a waveguide not only depensd on its intrinsic geometrical characteristics, but also depends on the waveguide geometry where it is inserted. 



§  CONCLUSIONS


In this paper we have presented a comprehensive method for extracting the effective polarizabilities of metamaterial elements patterned in 1D rectangular and 2D parallel plate waveguides. In our approach, the S-parameters—along with the knowledge of the normalized fields inside the waveguide—are used to find the electric and magnetic polarizabilities, making this technique a powerful tool for the design and characterization of metasurfaces. To demonstrate the validity of our approach, the polarizability extraction method was compared with direct extraction of the tangential components of the scattered fields, and excellent agreement between the two methods was demonstrated. We have also shown that the dipole modes predominate the scattering from metamaterial elements, since all higher-order multipole fields decay more rapidly with distance. The concepts presented in this paper pave the way for a simple and efficient approach to the analysis of metasurface antennas. The combination of polarizability extraction techniques with the dipole model provides a powerful and inherently multiscale modeling tool to design and characterize metasurface structures without any limitation on the element geometry or periodicity assumptions common to other homogenization techniques. 



§ RADIATION REACTION INSIDE A RECTANGULAR WAVEGUIDE 

 As described in section <ref>, the radiation reaction in a rectangular waveguide can be found by taking the real part of the surface integral of the Pointing's vector 𝐒. Let us consider again a thought experiment where two collocated electric and magnetic dipole 𝐩 and 𝐦 are placed at position 𝐫_0,in an environment that is described by the Green's function 𝐆(𝐫; 𝐫_0). The surface integral of the Poynting's vector is given by

    ∫𝐒· d𝐚 = iω (𝐩^*·𝐄 -  μ_0 𝐦^*·𝐇)

 where 𝐄 and 𝐇 correspond to the total fields inside the waveguide. Evaluating such fields at the dipole's location we get

 
    ∫𝐒· d𝐚=iω (𝐩^*·𝐄 -  μ_0 𝐦^*·𝐇) 
    
     = iω (𝐩*·𝐆_ee(r_0,r_0)·𝐩 -  μ_0 𝐦^*·𝐆_mm(r_0;r_0)·𝐦)


Now, considering the real part of Eq. <ref> it is possible to obtain a direct relationship between the total power radiated and the imaginary components of the Green's functions as

 
    Re{∫𝐒· d𝐚} = ω |𝐩|^2 Im{ G_ee(r_0,r_0) }
     
    - μ_0 ω |𝐦|^2 Im{ G_mm(r_0,r_0) }


On the other hand, from the modal expansion of the fields Eq. <ref> and Eq. <ref>, the same integral shown in Eq. <ref> results in

    ∫𝐒· d𝐚 =1/Z_n(|A_n^+|^2 + |A_n^-|^2).

where Z_n, A_n^+ and A_n^- have been previously defined in <ref>. By using the equations for the amplitude coefficients shown in Eq. <ref> into (Eq. <ref>) we obtain


    (|A_n^+|^2 + |A_n^-|^2) = ω^2 Z_n^2/4 (|E_n^+|^2 |p|^2 - μ_0|H_n^+|^2 |m|^2 )


In addition, it was previously demonstrated that for the fundamental mode, the fields E_n^+ and H_n^+ are normalized by means of Eq. <ref>. Replacing the specific expressions for the normalized fields, and equating Eq. <ref> and Eq. <ref> we obtain a direct relationship between the amplitude of the normalized modes and the imaginary Green's functions, which are our ultimate goal in this derivation. More explicitly,


    ω |𝐩|^2 Im{ G_ee(r_0,r_0) }  
    - μ_0 ω |𝐦|^2 Im{ G_mm(r_0,r_0) }
    
    = ω^2 Z_n^2/4 (|E_n^+|^2 |p|^2 - μ_0|H_n^+|^2 |m|^2 )


It can be observed that, in order to satisfy Eq. <ref> the terms multiplying the magnitude of the dipole moments |p|^2 and |m|^2 must be equal. After some algebraic derivation it is possible to conclude that

    Im{ G_ee(0) } = k^2/β ab  Im{ G_mm(0) }= k/ab.
 

