FDTD modeling of nonperiodic antenna located above metasurface using surface impedance boundary condition

This paper investigates an FDTD modeling method for precisely calculating the characteristics of a single, that is, a nonperiodic antenna located above a metasurface that consists of an infinite periodic conducting element on a flat dielectric substrate. The original FDTD method requires enormous computational resources to analyze such structures because an appropriate periodic boundary condition (PBC) is not supported, and a brute force approach has to be used for this reason. Another option is to use the array scanning method in which a single source is synthesized from a superposition of infinite phased array of point sources. In this method, some problems such as a mutual coupling between the single antenna and the metasurface, a computational error contained in a numerical integration over the Brillouin zone and so on have not been resolved yet. In order to resolve these difficulties and to reduce computational resources, a surface impedance boundary condition (SIBC) is incorporated into the FDTD method in this paper. The validity of the method is numerically confirmed by calculating an input impedance and a radiation pattern of a horizontal dipole antenna located above the metasurface.


Introduction
The use of periodic artificial materials or metamaterials has been investigated for antennas, microwave devices and many other applications in the field of optics due to their own high ability to enhance the original characteristics and to control electromagnetic wave propagation [1][2][3]. An example is a periodic leaky-wave antenna that a direction of antenna beam is changed by an exciting frequency keeping at an enhanced directivity. One of other applications is an artificial magnetic conductor (AMC) that can be used to create low profile antennas [4]. The AMC is often realized by a mushroom type electromagnetic band gap (EBG) structure that is used to suppress surfacewave propagation as well [5]. The purpose of this paper is to provide a computational method of electromagnetic fields excited by a single(nonperiodic) antenna that is placed near a two dimensional periodic structure. This type of antenna is sometimes referred to as a metamaterialinspired antenna in the field of antenna engineering.
The finite difference time domain (FDTD) method [6][7][8][9] has been widely recognized as one of the most powerful electromagnetic computational methods. Since the computational algorithm is simple and straightforward, and furthermore a practically acceptable accuracy can be easily obtained, the FDTD method has been applied to analyze the metamaterial antennas as well. For a perfectly periodic geometry such as a plane wave scattering by the metasurface, calculation costs can drastically be reduced using a periodic boundary condition (PBC), that is, Floquet's condition or Bloch's condition because the field computations are completed within a unit cell only [8][9][10]. On the contrary, the PBC can not be used when the single, that is, the nonperiodic antenna is included in the periodic structure. Therefore, a sufficiently large computational space that contains both quite a lot of periodic elements and the nonperiodic antenna is required to precisely analyze the electromagnetic fields when the original, that is, the brute force FDTD technique is used. Accordingly the computational resources increase remarkably in this brute-force FDTD(BF-FDTD) scheme. In order to overcome this difficulty, the array scanning method [11] has been incorporated into the FDTD calculation in which the single source is synthesized from a superposition of infinite phased array of point sources [12,13]. Although this method is theoretically attractive, an integral representation over the Brillouin zone are included in an electric field expression and its numerical calculation can not be accurately performed because an integrand exhibits a severe singularity. Furthermore, it may be extremely difficult to extend it to the method that can precisely calculate a mutual coupling between the single antenna and the metasurface. On the other hand, accurate surface impedance expressions for various metasurface have already been given, or can be easily calculated by the FDTD method or other numerical methods as well. Considering these situations, a computational scheme that the metasurface is replaced by the surface impedance and then a surface impedance boundary condition (SIBC) is incorporated into the FDTD calculation [14][15][16], is effectively available. In this paper, the validity and the accuracy of the SIBC-based FDTD method are numerically tested in calculations of an input impedance and a radiation pattern of a horizontal dipole antenna located above the metasurface, because this method is fundamentally an approximation method.
In the next section, a transmission line approximation of the surface impedance of the metasurface is briefly reviewed. These results are interpolated by a rational function of an angular frequency in order to successfully introduce it into the FDTD algorithm for analyzing a scattering problem due to a frequency dispersive medium. The metasurface treated here consists of the infinite periodic rectangular patch which is placed on the dielectric slab backed by the perfect conductor. The FDTD formulation of the electromagnetic fields for the frequency dispersive surface impedance is discussed in Section 3. In the FDTD method, a far-field is usually calculated using a far-zone transformation formula in time domain, but this paper uses the equivalence theorem in frequency domain directly. In Section 4, the computational scheme of the radiation pattern of the dipole antenna located above the metasurface is indicated after an exact electric field representation in far-zone is derived. It is numerically demonstrated using some numerical examples that the SIBC-FDTD is extremely effective for the far-field calculations.

Surface impedance boundary condition
An impedance boundary condition on a curved surface S illustrated in Figure 1 is expressed as [17,18] or simply E ¼ Z s ⋅K e [19], wheren is a unit normal to the boundary surface S, K e ¼n Â H and K m ¼ E Ân are an electric current and a magnetic current on S, respectively. Z s is a dyadic surface impedance and satisfies a condition n ⋅Z s ¼ Z s ⋅n.

Surface impedance of metasurface composed of rectangular patch on dielectric slab
The metasurface considered in this paper is illustrated in Figure 2. A W x Â W y PEC patch is periodically placed on the dielectric slab backed by the PEC. The periods in x and y directions are p x and p y , respectively. The surface impedance Z s of the flat structure as shown in Figure 2 can be represented by a scalar quantity in accordance with a polarization of an incident plane wave and has already been derived using the transmission line approximation [20], and is given by Z À1 s ðu ; vÞ ¼ Z À1 g ðu; vÞ þ Z À1 d ðu; vÞ where Z g (u, v) corresponds to the impedance for periodic patch and Z d (u , v) the impedance of a shorted transmission line of a length d. u is an incident angle as shown in Figure 2 and v is the angular frequency. When m r = 1, where When the approximate expression of the surface impedance is not explicitly given, the equivalent numerical data can be available by calculating the frequency and angular dependencies of the reflection coefficients for the plane wave incidence using the FDTD method or other computational techniques.

Rational function interpolation of surface impedance
A time domain surface impedance is needed in order to incorporate the SIBC into the FDTD computation, because the electric field is related to the magnetic field through the frequency dependent surface impedance as discussed in the above section. However, it is quite difficult to transform the surface impedance Z s (u, v) to the time domain impedance Z s (u, t) using equations (2)-(4), because Z s indicates a complicated frequency dependence.
On the other hand, the electromagnetic scattering investigated in this paper is considered as a kind of passive system such as a passive electric circuit. Therefore, it is reasonable that the impedance is expressed by a positivereal rational function with respect to the angular frequency v like a Foster function for a loss-less electrical circuit. In order to obtain the time domain surface impedance expression that is easily applicable to the FDTD method, we assume that the surface impedance is expressed as where s = jv, and A n and p n , (n = 1, 2, … , N) are unknown coefficients to be decided from the given data of Z s (s) for a fixed incident angle u [21]. The time domain surface impedance can be easily obtained from eq. (5), and is given by In order to confirm the validity and the accuracy of the interpolation method for the surface impedance given in Section 2.1, the coefficients A n and p n were decided by solving the least square equations obtained from given data of Z s (jv) from 10 GHz to 20 GHz. In this calculation, we 3 mm and e r = 4.0, and then 3 appropriate poles were found, that is N = 3 the approximated surface impedance and the original surface impedance are shown in Figure 3. It is found that the approximated surface impedance agrees with the original surface impedance perfectly.
3 Impedance calculation of dipole antenna using SIBC-FDTD

SIBC-FDTD formulation
In order to make full use of the capabilities of metasurface, the antenna is naturally placed very near the surface. Hence a height h shown in Figure 4 is small compared with an operating wavelength in practice. In this section, the FDTD method using the SIBC is studied for such situation. On the other hand, the radiated fields can be considered as the superposition from various angular spectrums and are incident to the metasurface with various angular spectrums. The reflect-backed field components to the antenna may be confined within the partial components incident to an area just under the antenna, so that the surface impedance at normal incident dominates the antenna impedance. Consequently, only the surface impedance at normal incident is incorporated into the FDTD formulation.
Transforming the SIBC E ðvÞ ¼ Z s ðv; 0Þẑ Â H ðvÞ to the time domain, next two convolutional representations are obtained at the boundary surface shown in Figure 4 E x ðtÞ ¼ À where Z s (t) is the Fourier transformed time domain surface impedance and is given by equation (6) in the proposed SIBC-FDTD calculation. Since the Z s (t) is expressed by an exponential function, the convolution integrals (7) can easily be evaluated by the recursive convolution (RC) method [22] or the piecewise linear RC (PLRC) method [23]. In this paper, the RC method was used for simplicity. The electromagnetic fields in other computational region are calculated by the conventional Yee's algorithm.

Numerical examples
In this section, the validity of the impedance calculation using the SIBC-FDTD is confirmed. The overall structure of the metasurface and its unit cell are the same as Figure 2, and the same parameters given in the above section were also used for all numerical examples. Above the interface at z = 0, the very thin horizontal dipole antenna whose length  is 2l = 10 mm, is placed at a position of z = h as shown in Figure 4. In all FDTD calculations shown below, the cell size was set Dx ¼ Dy ¼ Dz ¼ 0:1 mm. A PML (perfectly matched layer) absorbing boundary [6][7][8][9] was used for truncating the computational space. The SIBC-FDTD is compared with two different kinds of the brute-force FDTD methods. One of which is an infinite slab model shown in Figure 5a in which the dielectric slab is equivalent to infinite but the patches outside the PML are not included in the calculation. Other the finite slab model is shown in Figure 5b. This model is often used in the original FDTD and/or the commerciallyavailable simulators that are not provided with an absorbing boundary capable to a piecewise homogeneous medium as treated in this paper.
Since the computational space is truncated by the PML, we checked how large computational space is required for obtaining the accurate input impedance. The input impedance for the number of patches including the computational space is shown in Figure 6 where the total number of patches is given by N x Â N x and h ¼ 10 mm . 'BF ∞ ' denotes the brute-force FDTD using the infinite slab model shown in Figure 5a, and 'BF f ' is the one shown in Figure 5b. The operating frequency was set 15 GHz at which the surface impedance resonates as shown in Figure 3 and hence the metasurface acts as a perfect magnetic conductor (PMC) for a normal incident plane wave. It is found that 'BF ∞ ' converges approximately to the SIBC at N x = 30 which corresponds 3 wavelengths at 15GHz. On the other hand, the 'BF f ' slowly converges to the other value near 'BF ∞ ' in less than about 10 V. This relatively small difference may be caused by the reflected surface wave at the boundary between the dielectric slab and the air.
The frequency characteristics of the input impedance are shown in Figure 7 for h = 10 mm, 7.5 mm and 5 mm. In all BF-FDTD calculations, 30 Â 30 unit cells are included in the computational space. It is found that both methods agree fairly well with each other and that the impedance calculation by the SIBC method is effective even if the antenna is very close to the metasurface.

Radiation patterns 4.1 Far field representation
In order to investigate the scattering field from the metasurface, let us start with the exact field representation. The dyadic Green's functions for the dielectric slab backed by the PEC as shown in Figure 8 can be easily derived by the same way described in [24] or [25], and some formulas expressed by vector wave functions in cylindrical coordinate can be obtained. However, these expressions are inconvenient to use the scattering problem due to the rectangular patches shown in Figure 2. In this case, the rectangular coordinate representations of the dyadic Green's functions are available [26]. Using the Green's functions, the electric field in the region of z > 0 is expressed by   is different from the definition of the magnetic type Green's function in [24].
The far-field electric field is obtained by applying a saddle point method to equation (8), and is expressed as follows where D 0 (u, f) is a vector directivity function in free space and can be calculated in a familiar way. On the other hand, D s u; f ð Þ ¼ D wherek s ¼ sinu cosfx þ sinu sinfŷ À cosuẑ.
Although the far field is usually calculated using the farzone transformation formula in the time domain, a diffetent manner was adopted for efficient calculations in this paper. First, the antenna is excited by a sinusoidal voltage and then the amplitude and the phase of the equivalent currents K ðAÞ e , K ðAÞ m and J ðP Þ e are calculated at all points on the surfaces. The directivity function for an observation angle (u, f) is obtained by calculating the integrals (11) using the complex representations of equivalent currents.

SIBC approximation
When the SIBC is used, R (TE) and R (TM) are replaced by the next equations and the second term of D ðeÞ s in equation (11) is omitted.
where Z ðTEÞ 0 ¼ Z 0 =cosu, Z ðTMÞ ¼ Z 0 cosu, and the negative sign is adopted for the TM incidence.

Numerical examples
In order to confirm the validity of the SIBC-FDTD method for the far-field calculation, the radiation patterns for the same structure described above, were calculated using the BF ∞ -FDTD method and SIBC-FDTD method.
The radiation pattern D u (u) in z À x plane is shown in Figure 9 for the case that the length of the antenna 2l is half-wavelength at 10 GHz and the height h is quarterwavelength. A phase of the reflection coefficient is roughly estimated as 75 degrees from Figure 3, and hence the pattern is closely resemble the half-wavelength dipole antenna in shape even though h is is quarter-wavelength. The number of patches including in the computational space of the BF ∞ -FDTD method was set 60 Â 60 which is 4 times larger than the impedance calculations. Both patterns are normalized by the maximum value of the SIBC-FDTD. It is found that the SIBC-FDTD agrees completely with BF ∞ -FDTD, and that the validity of the SIBC is numerically demonstrated. Figure 10 shows the radiation pattern at 15GHz where 2l and h are half-wavelength and quarter-wavelength at this frequency respectively. In this geometry, the metasurface acts as the PMC and then a direct field from the antenna is reduced by the reflected wave from the metasurface near u = 0. It is found that the noticeable discrepancy between the BF ∞ -FDTD and the SIBC-FDTD is observed. In order to investigate the cause, the surface current density along x axis on the interface at z = 0 is shown in Figure 11. The discrete peaks correspond to the  current density on the conducting patch surfaces. Thus the current on the patch attenuates very slowly with respect to the distance from the antenna at the resonant frequency (15 GHz) but relatively rapidly decays at non-resonant frequency (10 GHz), and hence larger computational space is needed in the BF ∞ -FDTD calculation at 15GHz. Similar numerical error characteristics are hold even when the antennais extremely close to the interface as shown in Figure 12. It is often that the BF f -FDTD method is applied to the far-field calculation in expectation of getting practically acceptable results as is the same with the input impedance. In order to confirm it, the radiation pattern of the same antenna as Figure 10 is shown in Figure 13. It is found that BF f -FDTD method is obviously incorrect in spite of the huge computational space. This failure is caused by a large reflection from the edge of the dielectric slab. It is also suggested that the surface wave is effectively excited in the dielectric slab and has an extremely large effect to the farfield.

Conclusion
The SIBC-based FDTD method has been investigated to calculate the properties of the dipole antenna located above the metasurface that consists of the infinite periodic conducting element on the flat dielectric substrate. First, the rational function interpolation of the given surface impedance data was introduced in order to obtain the expression that can be successfully incorporated into the FDTD method. Next, the input impedance of the horizontal dipole antenna located above the metasurface which consists of the rectangular PEC patch element, was calculated using the SIBC-FDTD method and two kinds of BF-FDTD methods, that is BF ∞ -FDTD and BF f -FDTD. It has been shown that the SIBC-FDTD agrees fairly well with the BF ∞ -FDTD even if the computational space is relatively small and that the numerical accuracy of the BF f -FDTD method is interior to the BF ∞ -FDTD method but is acceptable. Finally, it has been shown that this method can be applied to the radiation pattern calculation successfully by enlarging the computational space. It has also been shown that BF f -FDTD is inappropriate because the the unecessary large reflection occurs in the dielectric slab.
Since the computational results of the BF-FDTD method depend on the number of the unit cell contained in the computational space, it must be checked for every metasurface structure how large computational space is required for obtaining the accurate results, especially the radiation pattern. On the contrary, the scattering fields by the surface impedance sheet may hardly be affected by the Fig. 11. Surface current density on interface along x axis.   dimension of the computational space even if the appropriate absorbing boundary such as the PML is used. This is an additional advantage of the SIBC-FDTD method. However, only the surface impedance at normal incidence is incorporated to the FDTD impedance calculation and the allocated positions of the magnetic field are different from the position of the electric field. Hence, the numerical error resulted from these inherent features are contained in the proposed SIBC-FDTD method. It would be the remaining problems to be solved for more accurate calculation. The paper has been treated the single dipole antenna, but is easily applicable to any other finite size antennas.