Comparison of Methods for Calculating the Field Excited by a Dipole Near a 2-D Periodic Material

A comparison of methods for calculating the field from a single dipole source in the proximity of an infinite periodic artificial material (PAM) is given. A direct plane-wave expansion method (PWM) is compared with the ldquoarray scanning methodrdquo (ASM). The ASM is shown to be the most efficient method for sources that are close to the PAM, since it only requires integration in the wavenumber plane over the Brillouin zone rather than over the entire wavenumber plane. It is also shown how the ASM may be used to efficiently calculate the Fourier transform of the field at any aperture plane, which is useful for asymptotic analysis. The ASM is shown to be more efficient than the PWM or the reciprocity method for calculating the transform of the aperture field. A discussion of the nature of the singularities in the complex wavenumber plane is given, and numerical issues associated with the implementation of the ASM are discussed.


I. INTRODUCTION
P ERIODIC artificial materials (PAM) have been the subject of great interest and have been increasingly applied to antennas and microwave devices [1]- [11].As an example, electromagnetic band gap (EBG) materials are used to suppress surface-wave propagation [3]- [8].In other applications, placing an antenna near/within a PAM has been shown to significantly affect the input impedance or to create directive beams [9]- [11].Additionally, many leaky-wave antennas are in fact periodic in Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.
The present contribution is focused on the efficient calculation of the field due to a single (non-periodic) dipole source in the presence of an infinite three-dimensional PAM that is periodic in two dimensions and finite in the third.Some of the related numerical implementation issues are also discussed.The PAM under consideration is allowed to have an arbitrary metallization within the unit cell.Although the general case of a skewed lattice could be treated, the discussion is limited here to a rectangular lattice with periods and along the and axes, respectively.A typical periodic structure (with a metallic cross element) is shown in Fig. 1, which shows the dipole source location and the observation point .The methods discussed can easily be applied to layered structures having dielectric layers, but the discussion here will focus on structures in free space for simplicity.
It is first reviewed how the field produced by a dipole source in the presence of the infinite PAM may be calculated by a direct plane-wave expansion method (PWM).Although this method is well known and (as shown here) is not always an efficient method, it is useful as a benchmark for comparison.
The "array scanning method" (ASM) for calculating the field due to a single dipole source near an infinite PAM is then reviewed.This method has been introduced previously and used in the analysis of phased arrays [14]- [16].The method has also been discussed recently in the context of analyzing the fields near one-dimensional (1-D) periodic artificial material structures [17].The purpose here is to extend the investigation of the ASM to two dimensions, and to provide a comprehensive comparison of the method with other methods for analyzing two-dimensional artificial materials.The focus is on providing a general methodology that allows the ASM to be incorporated into general moment-method integral equation solvers.This paper represents the extension to the 3-D case (periodic in two directions and finite in the third direction) of the 2-D case (periodic in one direction, infinite in another direction, and finite in the third direction) treated in [18] and the same formalism is used here.
In addition to providing a very efficient method for calculating the field of a dipole source, it is shown how the ASM may be used to efficiently calculate the Fourier transform of the field at any aperture plane of interest.This is useful for performing asymptotic analysis, since such analyses often start with the assumption that the relevant integrals that describe the field are in the form of an inverse transform integral.Having the field expressed in this form is also useful for complex path deformations to identify the types of wave phenomena that are present and to determine the launching amplitude of surface and leaky waves.It is shown that the ASM provides a more efficient calculation of the transform of the field than does the PWM or an alternative reciprocity-based method that is introduced.
In the ASM method the PAM is excited by an infinite periodic array of dipole sources with the same period as the PAM and phased with variable wavenumbers .An integration in (i.e., scanning the phased dipole array) over the Brillouin region, and , recovers the field produced by the single dipole source near the PAM.
An understanding of the nature of the and planes in the array scanning method is important for an asymptotic evaluation of the field along the interface of the PAM.It is shown that due to the periodic nature of the problem, there is an infinite set of branch point loci.Furthermore, there is also an infinite set of periodic pole loci that lead to surface and leaky waves.A discussion of these issues is given so that numerical aspects in the implementation of the ASM can be understood.
In Section II the PWM is reviewed.In Section III the array scanning method is discussed.In Section IV a comparison is given between the ASM and the direct PWM, and it is shown how the ASM field representation may be transformed into the PWM form.A comparison of computational efficiency between the PWM and the ASM is also given.The ASM representation is much more efficient when the source is close to the PAM, and this is discussed.In Section V it is shown how the solution from the ASM may be used to cast the solution in the form of a 2-D inverse Fourier transform integral, allowing for the efficient determination of the Fourier transform of the aperture field.In Section VI issues associated with the numerical implementation of the ASM are discussed, including singularities in the complex wavenumber plane and errors due to numerical integration.In Section VII numerical results are given that illustrate some of the features discussed in the previous section.Finally, in Section VIII conclusions are given.

II. SYNTHESIS BY PLANE-WAVE SUPERPOSITION
Consider the problem illustrated in Fig. 1 where a PAM is excited by an impressed unit electric dipole current source at , oriented along the direction (which may be or any other direction), described by where .Here and in the following, boldface symbols are used to identify vector quantities and the caret " " identifies unit vectors.The PAM structure in Fig. 1 is periodic along and , with periods and , respectively.The structure may also be periodic along , truncated after a finite number of layers, thereby constituting a slab of artificial material.In this case the term "supercell" is used to denote a unit cell in the and directions (which contains multiple cells spaced along ).For simplicity we consider here only metallic scatterers such as those shown in Fig. 1.In what follows, we represent a general observation point (which may be located in the -th unit cell) as , where is assumed to lie within the (0,0)-th unit cell.The field in free space at radiated by the single dipole source at is denoted by and it can be represented in terms of plane waves as where , and is the transverse part of . is the identity dyad.The "minus (plus)" sign is used when the observation point is above (below) the source point.Assuming the dipole source is above the periodic structure for simplicity, the incident field from the dipole source that illuminates the periodic structure is a continuum of plane waves of the form (4) with (5) Each incident plane wave in the spectrum is polarized in a direction where (6) (the asterisk denotes complex conjugate).The scattered field due to a unit-amplitude incident plane wave with transverse wavevector is denoted as .(A bar over a quantity is used here to signify that the quantity is either a unit-amplitude incident plane wave or is produced by such an incident wave.)Similarly, the total field due to a unit-amplitude incident plane wave is denoted as .From superposition, the scattered field produced by the dipole source in the PAM environment is evaluated by summing all the fields produced by each plane wave in (2).The scattered field produced by the single point source is thus (7) and similarly for .The electric field produced by a unit-amplitude incident plane wave is found by first solving the electric field integral equation (EFIE) within the supercell to find the surface current on the elements within the supercell.The field scattered by the PAM is then given as (8) The free-space periodic Green's function is (9) where the dyad is given by (3) with and the wavevector of each Floquet wave is , with transverse wavevector (10) where (11) (12) The longitudinal Floquet wavenumber along is (13) where is the free-space wavenumber.As will be clear later, it is convenient to exhibit the singularity explicitly in (9).Branch points are defined from , corresponding to .Therefore, if is fixed, there is a doubly infinite set of branch points in the plane at (14) If the integral is performed first in the double integration of (7), as indicated by the integration order, then the path of integration lies in the complex plane for each value of .The path stays on the top Riemann sheet of all the branch points, and the wavenumber is chosen as either purely real and positive or purely imaginary with a negative imaginary part (assuming a lossless medium).Geometrically, the location of the branch points for a given value of is determined by the intersection of a horizontal line (constant ) with the circles shown in Fig. 2.
The field radiated by the scattering structure when illuminated by a unit-amplitude incident plane wave can be expressed as a Floquet expansion, so that (15) Insertion of ( 9) into ( 8) yields ( 16) A significant numerical effort may be required to evaluate the scattered field from this equation since it involves a 2-D infinite integral and a 2-D infinite sum.This is particularly true when the source point lies vertically close (in the direction) to the PAM, since the weight function then contributes minimally to the convergence of the integrations for large values of .In fact, as the source dipole approaches the PAM, the weight function defined in (5) does not tend to zero as .However, in a periodic moment-method solution, the integral of the incident electric field over the basis functions on the unit supercell would be required for the computation of the right-hand elements, and this integration (which results in the Fourier transform of the basis functions) would result in a converging integral expression.
For source points that are well removed from the PAM, the integration converges rapidly, and is limited to essentially the visible space region .The ASM for calculating the scattered field, which remains efficient even for sources close to the PAM, is discussed in the next section.

III. THE ARRAY SCANNING METHOD (ASM)
The single dipole source in (1) oriented along the direction can be realized by using the "array scanning method," [14]- [16], i.e., by synthesizing the single dipole source from an infinite phased array of identical point sources located at , as shown in Fig. 3. Mathematically, the phased array of dipole sources is represented as (18) The wavenumbers and are the phasing gradients along the and directions, respectively.According to the ASM [15], the current of the single dipole can be obtained by integrating (18) over the Brillouin zone shown in Fig. 2, so that (19) Physically, this represents the fact that when the phased array currents are integrated in over the Brillouin zone, all of the currents located at in the phased array integrate to zero except the one located at .The incident field generated by is represented as , where the periodic dyadic Green's function is given in (9).The incident field is thus rewritten as (20) with (21) The incident field exhibits an infinite set of -indexed branch loci in the plane defined by the equation for that, for fixed , define the doubly infinite set of branch points in (14).The scattered field can then be represented as a sum of scattered Floquet waves of the form (22) where ( 23) and is the current on the scatterers in the (0,0) supercell, when excited by the infinite phased array of dipoles.This current is obtained from a numerical solution of the EFIE on the (0,0) supercell using a periodic moment-method code.The current is a periodic function of the spectral variable , with periods and in and , respectively, since the phased array of dipole sources is periodic in and .
According to the ASM, the current on the PAM is obtained from the infinite phased-array current by a spectral integration over the Brillouin zone, so that (24) The scattered and total fields produced by the single dipole source are also obtained by a spectral integration over the Brillouin zone, so that The integrands and are both periodic functions of , and with periods and , respectively.The singularities of the integrands of (25) and ( 26) are discussed in Section VI, and they are important to determine the proper numerical treatment of the integrals (see [17] for a similar problem with 1-D periodicity).
One important property of the ASM is that it only involves the calculation of a field from a periodic source (for either the incident or scattered field), and for this type of calculation highly efficient acceleration schemes may be used, such as the Ewald method [19].In fact the Ewald method may be used to calculate both the incident and scattered fields, instead of using ( 21) and (23).The Ewald method represents the periodic free-space Green's function as a sum of two double series, each of each has Gaussian convergence.Typically, each double series in the Ewald method converges with summation limits from to 2. In all of the numerical calculations performed later, the Ewald method was used.
In the expressions for the surface current and the electric fields and we can assume for convenience and without loss of generality that the point is located within the (0,0) unit cell and that the observation point is now described by .In this case expressions (24), (25), and (26) can be rewritten as (29) where denotes either , or .

A. Equivalence of Methods
By noticing that the incident field in (20) produced by the infinite phased array of dipoles in the ASM can be viewed as a weighted superposition of plane waves, the relation between the scattered field produced by and the field produced by a unit-amplitude incident plane wave may be written as The scattered field calculated by the ASM may then be cast into the form (33) Mathematically, ( 7) is seen to be equivalent to (33) by using the fact that (34) in order to transform the wavenumber integrations from the entire real axis into one over the Brillouin zone.The integrand in (25) can in principle be calculated using (30).However, doing so requires the solution of an infinite number of plane-wave scattering problems.A much more efficient method is to directly calculate by numerically solving the problem of an infinite set of phased dipole sources above the PAM.This numerical solution requires a periodic moment-method analysis that is essentially no more numerically intensive than that for a single plane wave scattering problem.
The advantage of the ASM is that (25) requires a spectral integration that is carried out only over the Brillouin zone in contrast to (17), which requires an integration over the entire plane.(However, as mentioned previously, for source points that are well removed from the PAM, the integration in (17) will converge very rapidly.)

B. Numerical Comparison of Methods
To clarify the comparison of computational efficiency between the ASM and the PWM, both methods may be compared in terms of how many periodic moment-method solutions they require in order to determine the current on the PAM due to the original source dipole (the determination of the current on the structure is often the most time-consuming step).The details of this comparison are provided in the Appendix.A ratio is defined as the number of periodic moment-method solutions required in the PWM relative to the number required in the ASM.The ratio depends dramatically on the vertical separation , where is the coordinate of the PAM (assuming for the sake of simplicity that it is planar here).From the analysis in the Appendix, it is seen that for very large vertical separations, is roughly on the order of (35) However, as decreases the ratio increases, with the rate determined by the size of the testing basis/testing functions used in the moment-method solution for the PAM.
Table I gives the value of calculated from the rough orderof-estimate formula (55) for a typical case having .In this comparison the half-length of each testing function is assumed to be one of two values: either (representing a typical sub-domain basis function) or (representing a typical full-domain basis function).It is seen that for large separations the PWM has a slight advantage, while for small separations the ASM has a large advantage.According to this (rather crude) estimate, the ASM becomes more efficient for a vertical spacing that is less than about two wavelengths.For large vertical separations the size of the basis functions is inconsequential, but for larger separations it is obvious that using smaller basis functions results in a much slower convergence in the PWM.

V. FOURIER TRANSFORM OF APERTURE FIELD
The calculation of the transform of the field on the aperture plane of the periodic structure (or any other plane located at a fixed value of ) is important for performing asymptotic analysis, and for identifying the launching amplitude of surface and leaky waves.Therefore, this calculation is considered here.Three methods for calculating the transform are discussed and compared.

A. Plane-Wave Expansion Method
Starting from (17), the scattered field may be put in terms of an inverse Fourier transform integral by "collecting the spectrum."That is, all of the plane waves that have a wavevector are summed together.Mathematically, this is equivalent to using the change of variables and for each term of the sum in ( 17), followed by the substitution and , and results in the expression The scattered field thus has the inverse Fourier transform representation Considering the definition of the function in (5), expression (38) explicitly shows that the integrand has an infinite set of branch point loci, located at [see the discussion after (13)].
Although (38) is a valid expression for the Fourier transform of the field at any value, it may be numerically intensive to compute, since it requires the numerical solution of an infinite number of scattering problems involving an incident plane wave (i.e., an infinite number of incidence angles corresponding to wavevectors .)However, as mentioned previously, the weight function decays rapidly with increasing and when the source point is well separated from the PAM, so that expression (38) eventually becomes quite efficient for large separations.In the next section an efficient method for obtaining the Fourier transform using the ASM is presented, which remains efficient even for small separations between the PAM and the source point.

B. Array Scanning Method
It is possible to "unfold" the integration path from the Brillouin zone shown in Fig. 2 to the entire plane.Doing so allows for a convenient identification and calculation of the Fourier transform of the field at any fixed height .The path unfolding is done by first inserting ( 22) into (25), and then using the property that (39) The shift of variables and is applied for every term of the sum, leading to where is calculated using (23) and the supercell current . By comparing (40) with (37) the Fourier transform of the scattered field is identified as (41) Equation ( 41) indicates that the Fourier transform of the aperture field is (to within the simple constant ) the same as the amplitude of the (0,0) Floquet wave that is radiated by the PAM when it is excited by the phased array of dipole sources.It is straightforward to extract the amplitude of this fundamental Floquet wave from the periodic moment-method solution with the phased array of dipole sources.
In the case of field evaluation at a point larger than all the -points of the scatterer (i.e., when for all ) a simple explicit expression for the amplitude of the scattered fundamental Floquet wave (and hence the transform of the scattered field) can be derived (an analogous result applies when is smaller than all the -points of the scatterer.)Using ( 40

C. Reciprocity Method
The reciprocity method may be used to provide an alternative calculation of the Fourier transform of the total field at any fixed plane.Assume that the transform of the component of the electric field ( or any other direction) is to be calculated at for transform wavevectors .An infinite "testing current" sheet of the form ( 44) is placed at .Dot multiplying the scattered electric field from the PAM due to the dipole source with the testing current and integrating over the plane at yields the Fourier transform of the component of the electric field at .By reciprocity [20], this is equivalent to calculating the electric field radiated by the testing current sheet (in the presence of the PAM) and sampling the component of the field at the source dipole location .The incident field radiated by the current sheet in ( 44) is (45) This field is a plane-wave field that impinges on the PAM, and the total resulting field is then sampled in the direction at to yield .The reciprocity calculation requires a plane-wave excitation of the PAM for each wavevector .In comparison, the ASM formula (41) requires the excitation of the PAM by a phased array of dipole sources.Both problems require the same level of computational effort.However, in the ASM, a calculation of the phased-array dipole problem with wavenumbers ranging over the Brillouin zone is sufficient to reconstruct the Fourier transform for all wavenumbers in the entire plane.This is because of property (39).
Of course, in the numerical calculation of the Floquet waves (the left-hand side of (39)), enough Floquet terms must be calculated to ensure that the wavevector extends over the desired range for which the transform is being computed.In view of ( 39) and (41), the Fourier transform is known over the entire plane once the phased-array scattering problem has been solved for all wavenumbers within the Brillouin zone.Hence, the ASM is much more efficient numerically than the reciprocity method in calculating the complete Fourier transform for all wavenumbers .If the transform is only desired for a limited range of wavenumbers, the reciprocity method is of comparable efficiency.

A. Singularities
In performing the integration in (25) or (26), branch-point singularities may be encountered on the real axis of the plane.There is a doubly infinite set of branch points at in the plane, given by ( 14) and shown in Fig. 4. The branch points are periodic along the real axis, with period .In the central Brillouin zone, defined as , there are an infinite number of branch points along the imaginary axis.There may be a pair of branch points on the real axis, depending on the value of the integration variable .When and there will not be any branch points on the real axis.(This means that in Fig. 2 the branchpoint circles do not intersect vertically, and the horizontal line passes through one of the vertical gaps between the circles.)Otherwise, there will be at least one pair of branch points on the real axis.
In many practical cases the frequency is low enough so that and , and there will then be at most one pair of branch points on the real axis within the integration region , and all the others will lie on the imaginary Fig. 4. Spectral k plane, for a fixed k value.The poles (the x points below the branch cut) and branch points are periodic in the k plane, with period 2=a.(a) case with ka < , or equivalently a= < 1=2, (b) case with ka > .axis.This is the case shown in Fig. 4(a).The integration in the complex plane in the region can be deformed off of the real axis to avoid the branch points on the real axis.
When , branch cuts corresponding to adjacent values of in ( 14) may overlap on the real axis as shown in Fig. 4(b).This will happen when .In this situation it is not possible to deform the path to avoid the branch points in the integration region when staying on the top Riemann sheet.Hence, this case requires careful integration along the real axis, especially in the region where the branch cuts overlap.
It is also interesting to note that the integrand of (25) is more singular at the branch points than is the integrand of (26).This is because of the physical fact that the scattered field from the source decays with distance more slowly along the PAM than does the total field.The decay of the scattered field is , which matches that of the incident field (the two must cancel on the conducting surfaces of the PAM).However the total field has a more rapid decay of far away from the source, which is typical of the behavior observed in any interface problem, including, e.g., a source over a dielectric layer.Hence, especially when evaluating the scattered field, care must be taken to integrate accurately near the branch points, or to deform around them if possible (for the case of Fig. 4(a).The behavior of the fields is verified numerically in the next section.
There may also appear surface-wave poles on the real axis if the PAM allows for the guidance of such waves.In the case of Fig. 4(a), where , deforming the path off the real axis will allow for the avoidance of such poles.In the case of Fig. 4(b), where , it is not possible to deform the path entirely off the real axis (staying on the top Riemann sheet).However, in this case all guiding modes must be leaky modes [21], because at least one Floquet wave is a radiating fast wave for all values within the Brillouin zone.Hence, the leakywave poles will lie off the real axis.Leaky modes may also exist in addition to surface waves, in the case of Fig. 4(a).(Leaky-mode poles are denoted with crosses in the figure .)The leaky modes may cause numerical difficulty if they are close to the real axis, but this is an unusual situation unless the PAM has been designed to form a leakywave antenna type of structure.
The above discussion of singularities in the complex plane is common to both the ASM and the PWM.In the PWM the integration is over the entire plane, and hence the integration in is over the entire axis instead of only the region in Fig. 4.

B. Numerical Integration Error
If the integrals (25) and ( 26) are approximated by sampling uniformly the integrand in the Brillouin zone, a physical interpretation of the numerical error can be given.Besides providing a simple physical interpretation of the error, the use of the "midpoint rectangle rule" of integration is chosen because it is especially accurate for periodic functions, and in particular for those cases where the integrand is smooth [22].Sampling the integrand at points (46) (47) with and , i.e., at the center of each of the intervals (i.e., the midpoint rectangle rule), leads to an interesting physical interpretation of the error.In particular, the error in approximating the field at a location due to a dipole source at is equivalent to summing the field produced by an infinite number of dipole source "images" located at for (the proof is omitted here).The integration approximation is thus satisfactory when the nearest images (those with ) are sufficiently far away, i.e., and are large enough to result in a large spatial field decay from the nearest images.One interesting ramification of this is the fact that if the medium is lossy, the error eventually decreases exponentially with an increasing number of integration points.This surprising result is consistent with the discussion in [22], which shows how the error in the simple midpoint rectangle rule often has exponentially small errors when integrating smooth periodic functions.The loss makes the integrand smooth over the integration region, since branch points are now located off the real axis.For lossless cases, adaptive integration schemes are possible as a way to improve the numerical integration, with a coarser sampling away from the singularities in the complex plane and a finer sampling close to the singularities.

VII. RESULTS
Results are obtained using a uniform sampling (rectangle integration rule) to perform the numerical integrations, as discussed in Section VI-B.The array of scatterers is made of resonant strip-dipoles with length , where is the freespace wavelength, and width , oriented along .The metal dipole in the (0,0) unit cell is divided into 10 rectangular subdomains of length and 9 rectangular rooftop basis functions are used.The array is a square lattice with element spacings .The source is an elementary -directed electric dipole with a unit amplitude ( Am) located at (the origin is at the center of the (0,0) metal strip dipole).The current (in amps) is calculated on the metal dipoles at a distance from the left end of each dipole.Results are presented for the current on dipole centered at , where is varied in the results below.In Fig. 5 we show the percent error in the current defined as , where represents the current on the metal dipoles obtained using many integration points , but keeping the number of basis functions fixed at nine.The error thus calculated is therefore a measure of the error in the numerical integration over the Brillouin zone.In Fig. 5(a) the three curves represent the percent error at three different supercell locations, corresponding to along the -direction.The current is evaluated via the ASM described in Section III, using (29) to evaluate the surface current on the metal dipoles.The ASM is performed by integrating over the real and axes, with the same number of equidistant sampling points in ( 46) and (47), corresponding to a total number of sampling points in the plane.Fig. 5 shows the above mentioned error versus .Because the Fig. 6.The same result as for Fig. 5, but for a lossy air medium that has a loss tangent of 0.1.periodic structure is in free space, the periodic Green's function employed in the moment-method solution for the currents has been evaluated using the Ewald method [19], which is much more efficient than a direct evaluation of (9).The Ewald method has been used to calculate the incident field from the phased array of sources in the ASM, and also to calculate the field from the periodic currents on the metal dipoles in the ASM, which is used to generate the matrix elements in the moment-method system of equations.
Fig. 5(a) shows that the relative error is larger for those metal dipoles that are further away from the source because the field produced by the source is weaker there and the effect of the nearest "images" at , for , is thus relatively stronger (see the discussion in Section VI-B).In Fig. 5(b) the percent error is plotted for (the metal dipole in the (0,0) cell), using Gaussian quadrature with three different orders ( and ).In this comparison equals the order times the number of subintervals used in the integration.Using order (one sample point per subinterval of integration) is equivalent to the midpoint rectangle rule of integration, and hence the curve is identical to the curve in Fig. 5(a).The results show that the error is about the same when using the midpoint rule as with Gaussian quadrature (see the discussion in Section VI-B for an explanation of why the midpoint rectangle rule is expected to work well).
Fig. 6 shows the same result as in Fig. 5, but for a lossy case where the air medium now has a loss tangent of 0.1.Comparing Fig. 7.The current on the metal dipoles calculated by using the ASM versus the number of cells away from the (0,0) cell in the y direction.(a) A result is shown using a large number of integration points.Also shown is a normalized 1=n curve for comparison.(b) Results are shown for various numbers of integration points.
Figs. 6(a) and 5(a), it is seen that the loss significantly increases the convergence rate of the numerical integration.This is because the lossless case of Fig. 5(a) had branch points on the real axis within the range of integration.Adding loss moves the branch points off the real axis, resulting in a smoother integrand.Equivalently, the addition of loss renders the field produced by the "images," the reason for the error when using the midpoint rectangle integration rule, weaker.The midpoint rectangle rule now outperforms Gaussian quadrature, as Fig. 6(b) shows.The error also decreases exponentially with the number of integration points, as expected from the discussion in Section VI-B.
In Fig. 7(a) the current is plotted over the first 100 metal dipoles along the direction ( varies from 0 to 100).The field decay follows the expected behavior as seen by the addition of the normalized curve.The curve has noticeable oscillations, which are not due to numerical error is sufficient to guarantee accurate results).Evidently, this is due to interference between the main space wave and a weaker mode (perhaps a leaky mode) propagating along the PAM interface.
In Fig. 7(b) the field is plotted again along the direction using an accurate solution that is compared with that obtained using and .Using fewer integration points causes loss of accuracy away from the source.For points the effects of the nearest "image" source in cell is clearly visible.The null at is caused by the cancellation of the field produced by the source and this image.

VIII. CONCLUSION
A comparison of methods for calculating the field due to a single dipole source in the presence of an infinite periodic artificial material (PAM) has been provided.A direct plane-wave expansion method (PWM) has been compared with the "array scanning method" (ASM), and the latter was shown to be much more numerically efficient when the source is located close to the PAM, since it only requires an integration over the Brillouin zone rather than over the entire wavenumber plane.
For some applications, such as asymptotic analysis and calculation of surface-wave and leaky-wave modal amplitudes, it is also important to be able to calculate the Fourier transform of the field at some aperture plane, so that the field is in the form of an inverse transform integral.Three methods were compared for the calculation of the transform: (1) the PWM, (2) the ASM, and (3) a reciprocity-based method.Of these, the array scanning method is the most efficient.The PWM requires the solution of an infinite number of scattering problems (i.e., different incident waves impinging on the PAM) in order to determine the transform for a given set of wavenumbers .The ASM and the reciprocity method require only a single scattering problem.However, the array scanning method allows for the complete transform to be constructed from the scattering solution for incident wavenumbers that are within the Brillouin zone.
The nature of the complex plane in the wavenumber integration of the ASM was discussed in some detail, including the location of the branch points and well as surface-wave and leaky-wave poles.An understanding of the singularities in the complex plane is important for numerical implementation.A physical discussion of the error involved in approximating the integrations with discrete samples was also given to help understand the numerical error that arises in the method.

Manuscript received June 9 ,
2006; revised December 29, 2006.This work was supported in part by the EU-funded project METAMORPHOSE (FP6/NMP3-CT-2004-500252) and in part by the State of Texas Advanced Technology Program.F. Capolino is with the Department of Information Engineering, University of Siena, Siena 53100, Italy (e-mail: capolino@dil.unisi.it).D. R. Jackson and D. R. Wilton are with the Department of Electrical and Computer Engineering, University of Houston, TX 77204-4005 USA.L. B. Felsen, (deceased), was with the Department of Aerospace and Mechanical Engineering and the Department of Electrical and Computer Engineering, Boston University, Boston, MA 02215 USA, and also with Polytechnic University, Brooklyn, NY 11201 USA.

Fig. 1 .
Fig. 1.Periodic structure made of a rectangular lattice of scatterers with dipole source at r and observation point at r.

Fig. 2 .
Fig. 2. Periodic branch-point loci in the plane Rek ; Rek .The area contained within the dashed line represents the Brillouin zone that is used in the ASM.

Fig. 3 .
Fig. 3. Periodic structure excited by an infinite set of phased dipole sources at locations r = r +max+nbŷ.The original dipole source is located at r .
be conveniently represented as a sum of Floquet waves, as (27) where (28) three-dimensional Fourier transform of the current on the unit cell.

Fig. 5 .
Fig. 5.The percent error in the current on the metal dipoles calculated by using the ASM versus the number of sample points along k or k , with the same number P = Q of equidistant sampling points in (46) and (47).(a) Error shown for three different (0,n) metal dipoles when using the midpoint rectangle rule of integration.(b) Error shown for the (0,0) metal dipole for three different orders of Gaussian quadrature (k = 1 is equivalent to using the midpoint rectangle rule).

TABLE I THE
RATIO R, GIVING THE COMPUTATIONAL EFFORT OF THE PWM RELATIVE TO THAT OF THE ASM, FOR VARIOUS VERTICAL SEPARATIONS OF THE SOURCE FROM THE PERIODIC STRUCTURE