IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 54, NO. 4, APRIL 2006 1619
Periodic FDTD Analysis of Leaky-Wave
Structures and Applications to the Analysis of
Negative-Refractive-Index Leaky-Wave Antennas
Titos Kokkinos, Student Member, IEEE,
Costas D. Sarris, Member, IEEE
, and
George V. Eleftheriades, Fellow, IEEE
Abstract—The determination of the attenuation constants of
periodic leaky-wave structures via the finite-difference time-do-
main (FDTD) method has been pursued so far via the simulation
of a number of unit cells that is large enough to guarantee the
convergence of the computed value. On the other hand, Brillouin
diagrams of periodic structures can be readily extracted via the
simulation of a single unit cell, terminated in periodic boundary
conditions. This paper demonstrates a methodology that enables
the concurrent extraction of leaky-wave attenuation constants and
Brillouin diagrams of periodic structures, through the FDTD sim-
ulation of a unit cell. The proposed methodology is first validated
and then employed to model leaky-wave radiation from a two-di-
mensional negative-refractive-index transmission-line (NRI-TL)
medium. Apart from evaluating the characteristics of forward and
backward leaky-wave radiation from such a medium, a lumped-el-
ement macro-model, with element values determined from the
FDTD simulation, is extracted. The FDTD analysis, combined
with this equivalent circuit, is used to investigate theoretically the
possibility of the NRI-TL medium, as a leaky-wave antenna, to
achieve continuous scanning from backward to forward end-fire
and broadside radiation.
Index Terms—Complex propagation constants, finite-difference
time-domain (FDTD) methods, leaky-wave antennas (LWAs), neg-
ative-refractive-index (NRI), periodic boundary conditions.
I. INTRODUCTION
L
EAKY-WAVE antennas (LWAs) have been extensively
studied for several years, because of their comparative
advantages, namely high directivity and frequency beam-scan-
ning, but also the richness of the electromagnetic phenomena
related to them. Recently, the study of leaky-waves in negative-
refractive-index (NRI) metamaterial-based structures has in-
dicated the possibility of designing novel types of one- and
two-dimensional LWAs, thus adding a new direction of research
to this area.
Theoretically studied by Veselago in [1], NRI media con-
currently exhibit negative dielectric permittivity
and
Manuscript received October 30, 2005; revised December 23, 2005. This
work was supported by the Natural Sciences and Engineering Research Council
of Canada under a Discovery Grant and a Strategic Grant.
T. Kokkinos was with the Edward S. Rogers Sr. Department of Electrical and
Computer Engineering, University of Toronto, Toronto, ON, Canada M5S 3G4.
He is now with the Wireless Communications Research Group, Loughborough
University, Leicestershire LE11 3TU, U.K.
C. D. Sarris and G. V. Eleftheriades are with the Edward S. Rogers Sr. Depart-
ment of Electrical and Computer Engineering, University of Toronto, Toronto,
ON, Canada M5S 3G4 (e-mail: cds@waves.toronto.edu).
Digital Object Identifier 10.1109/TMTT.2006.871367
magnetic permeability
and, as a result, a negative re-
fractive index
. The latter implies the support of back-
ward-waves (characterized by antiparallel phase and group ve-
locities or the formation of a left-handed triplet by the wave
vector
, the electric field vector , and the magnetic field )
by these media. Recent realizations of NRI “metamaterials.”
in the form of composite periodic structures of negative effec-
tive parameters, have enabled the experimental investigation of
their unconventional properties and salient features. While the
split-ring resonator and thin-wire structure of [2] offered the first
demonstration of an NRI medium, a loaded transmission-line
(TL) approach provided a planar and broad-band alternative im-
plementation [3]–[6].
One of the several loaded NRI-TL metamaterial-based appli-
cations that emerged in the literature is their use for the design of
one-dimensional (1-D) LWAs radiating in the backward end-fire
direction [7], [8] or having a beam-scanning capability from
backward to forward end-fire [9], [10]. Furthermore, two-di-
mensional (2-D) LWAs, based on loaded NRI-TL grids, have
been studied in [11]–[13]. Inspection of the dispersion diagram
of the 2-D loaded NRI-TL metamaterial grid of [4] (shown in
Fig. 1) reveals the existence of a fast wave region (where the
phase constant
of the structure is smaller than the free space
propagation constant
). Within this region, first backward and
then forward waves are supported. Then, the phase-matching
condition at the interface between the dielectric substrate of the
metamaterial and air enforces the emergence of a leaky-wave
beam at an angle
(1)
implying the operation of the metamaterial grid as a 2-D LWA.
Furthermore, the transition of
from negative to positive values
suggests that the angle of emergence of the leaky waves can
theoretically vary from
to .
The aforementioned structure is the motivating application
for this paper’s research into the full-wave analysis of LWAs via
the finite-difference time-domain (FDTD) technique. A partic-
ular question of interest that stems from the survey of previous
work in this area is whether it is possible to extract the attenua-
tion constant of leaky-wave periodic structures from the FDTD
analysis of their unit cell. This question becomes intriguing by
the fact that the translation of Floquet’s boundary conditions
from the frequency domain (where they are naturally cast) to
0018-9480/$20.00 © 2006 IEEE
1620 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 54, NO. 4, APRIL 2006
Fig. 1. Brillouin diagram of the loaded NRI-TL unit cell of [22], extracted
through the periodic FDTD analysis of [21], around the
=0
-region. The
results, for both the slow and fast waves, are in agreement with those extracted
through the TL theory-based analysis of [23] and Ansoft’s HFSS.
the time domain is invariably pursued by fixing the phase dif-
ference between the field components at the boundaries of a unit
cell to a real number and determining the resonant frequencies
of fields sampled within the cell [14]–[17]. On the contrary, pre-
vious computations of attenuation constants of LWAs were ex-
clusively made by simulating a relatively large number of unit
cells, terminated into regular absorbing boundary conditions,
until the value of the convergence of the attenuation constant
was achieved [29].
This paper builds on previous work by the authors [18], [19]
to demonstrate a methodology that allows the extraction of com-
plex propagation constants from an FDTD simulation of a unit
cell of a leaky-wave structure terminated in periodic boundary
conditions. The proposed technique is initially validated and
then applied to the modeling of leaky-wave radiation from the
loaded NRI-TL grid of [4]. In particular, the frequency depen-
dence of the attenuation constant of the forward and backward
fast waves in the structure is determined. This dependence is
considered when the stopband indicated in the dispersion dia-
gram of Fig. 1 at
is open and in the limiting case of
that being closed. Finally, the FDTD simulations allow for the
extraction of a lumped-element equivalent circuit that captures
the physical behavior of the structure and illuminates some of
its qualitative features.
II. L
EAKY-WAV E STRUCTURE MODELING IN PERIODIC
FDTD IMPLEMENTATIONS
This section is aimed at investigating the possibility to extract
complex propagation constants of leaky-wave periodic struc-
tures via the FDTD analysis of their unit cell. To that end, the im-
plementation of periodic boundary conditions (PBCs) in FDTD
is briefly revisited and carefully inspected.
Consider a 2-D periodic structure of periodicities
and
in the and directions, respectively. Then, its direct lattice
vector is
. Letting the reciprocal lattice wave
vector of the structure be
, Floquet’s theorem, in
the frequency domain, can be cast in the form of the following
equations:
(2)
(3)
where
and are phasors of the electric and magnetic fields,
respectively, while
is a
complex propagation constant, where
is the
corresponding attenuation constant vector.
In the past, the translation of (2) and (3) from the frequency to
the time domain and their incorporation into the FDTD update
equations was negotiated for the case of purely real propagation
constants [20]. In general, the direct implementation of these
equations in an FDTD framework is impossible, since it results
in update equations involving future (unknown) field values.
In order to circumvent this problem, several techniques have
been proposed. As mentioned in the Introduction, their common
starting point is to fix the wave-vector components to real values
and deduce the frequencies to which they correspond via the
computation of the resonant frequencies of the sampled fields. A
question that is naturally posed is whether such a methodology
would preclude the determination of any attenuation constants
involved, limiting the applicability of these techniques only to
cases where
.
In [21], the authors applied the sine–cosine method of [14] in
order to perform the dispersion analysis of the NRI medium of
[4]. The part of the dispersion diagram corresponding to axial
propagation (the
-part of the Brillouin zone) is shown
in Fig. 1. The PBCs used for this diagram are:
,
with
(in radians) and . Note that, al-
though the structure is 2-D, only axial propagation is considered
throughout this paper, allowing for the unambiguous use, hence-
forth, of
and . The diagram includes both slow
and fast waves, whose resonant frequencies were computed in
a unified way via the sine–cosine method. However, an inspec-
tion of the associated field resonances for slow and fast waves
indicates an important difference that is illustrated through the
following two examples. First, when the phase difference corre-
sponding to point A of Fig. 1, with
rad, is enforced
between the
-directed boundaries of a unit cell, then the ver-
tical electric field inside the structure exhibits the nondecaying
resonance pattern shown in Fig. 2. On the other hand, moving
the point of interest inside the fast-wave region of the structure,
when the phase difference corresponding to point B of Fig. 1,
with
rad, is enforced, the field pattern exhibits a
decaying resonance, as shown in Fig. 3. In both cases, the field
initially follows the source excitation (a Gabor pulse). When
the excitation fades away, the mode that corresponds to the en-
forced
is established within the cell. In the absence of any
ohmic or dielectric losses inside the unit cell, one can only inter-
pret this decaying resonance as evidence of radiation, success-
fully accounted for by the FDTD method, even though PBCs
without any explicit enforcement of an attenuation constant are
employed.
These numerical observations have an important implication
regarding the nature of methods that implement PBCs in FDTD,
KOKKINOS et al.: PERIODIC FDTD ANALYSIS OF LEAKY-WAVE STRUCTURES AND APPLICATIONS TO ANALYSIS OF NRI LWAs 1621
Fig. 2. Time-domain evolution of the vertical electrical field of a slow-wave
mode within the unit cell of the loaded NRI-TL metamaterial structure (termi-
nated with the PBCs
d
=
0
0
:
6
rad,
d
=0
rad).
Fig. 3. Time-domain evolution of the vertical electrical field of a fast-wave
mode within the unit cell of the loaded NRI-TL metamaterial structure (termi-
nated with the PBCs
d
=
0
0
:
07
rad,
d
=0
rad).
such as the split-field and the sine–cosine methods. When the
propagation constant
along the axis is fixed to a certain real
value, the condition that is numerically realized is NOT
since the latter would not allow for a decaying resonance, but
(4)
.
In order to support the validity of this statement, the field
update procedure in a periodic FDTD code is analyzed. For
simplicity, let us consider the one-dimensional (1-D) domain
of a unit cell (Fig. 4), which corresponds to a TEM-wave case.
Electric field nodes are positioned on the periodic boundaries.
Fig. 4. Updates of field components at the boundaries of the unit cell of a 1-D
periodic structure through the application of PBCs and FDTD update equations.
Fig. 5. Evolution of the electric field amplitude within a single unit cell of a
lossy structure (with PBCs enforced along the axis of periodicity) compared
with the one in three cells of the corresponding infinite periodic structure.
Then, the updates proceed as follows. Solving Maxwell’s equa-
tions, one can update
on node A, since its neighboring
is known from through the PBC (magnetic field updates
are complete, once the electric field updates start). Then,
is used to update , which is necessary for the magnetic field
updates that follow. Note that the direction of this update is
forward (a node at
is used to update a node at ). Once
the magnetic updates proceed,
is updated by Maxwell’s
equations, and its value is transferred to
and so on. Note that
the direction of this update is backward (a node at
is used to update a node at ). Therefore, if the modal
field that FDTD solves for is leaky, then the spatial field decay
cycle is not interrupted upon the enforcement of the PBCs, pre-
cisely because of this leap-frog sequence of boundary updates.
On the contrary, this cycle is extended and taken care of through
the implicit application of the boundary conditions in FDTD
and the solution of Maxwell’s equations, which produce a field
decay in space at a rate set by the attenuation coefficient. This
is the reason that (4) is the condition actually implemented in a
periodic FDTD technique. The way a leaky mode is accounted
for by a periodic FDTD technique is clearly depicted in Fig. 5,
where the spatial evolution of the field amplitude within a single
unit cell is shown, along with the field amplitude in three cells
of the infinite periodic structure, as a function of space. The key
observation is that the spatial attenuation of fields in an infinite
periodic structure is reproduced within the unit cell, during the
FDTD time-marching procedure, allowing for the calculation of
1622 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 54, NO. 4, APRIL 2006
Fig. 6. Field sampling points within a unit cell of a periodic structure for the
calculation of the complex propagation constant with (5).
its rate (which is the attenuation constant sought for). It is finally
noted that one would arrive at the same conclusion no matter
how the field nodes had been arranged in the domain of Fig. 4.
The outstanding question of
how the attenuation constant of
a leaky wave can be computed is addressed in the following
section.
III. M
ETHODOLOGY FOR THE
PERIODIC FDTD A
NALYSIS OF
LEAKY-WAV E STRUCTURES
A. Attenuation Constant Calculation
For the determination of the attenuation constant
in a pe-
riodic leaky-wave structure, a standard technique employed for
the calculation of attenuation constants in lossy/radiating (non-
periodic) guiding structures [25] is applied. However, the do-
main of application of this approach is now the unit cell of a
periodic structure and is terminated at PBCs that enforce, per
the analysis of the previous section, the phase constant of the
waves excited in the cell to that of a fast wave in the structure.
In particular, if
, which is the waveform of either an elec-
tric or a magnetic field component, is sampled at two points
and (Fig. 6) along the axis of periodicity, then the complex
propagation constant
can be deduced as
(5)
where
denotes a Fourier transform. Evidently, the real part of
the complex propagation constant
in (5) has to be the same
as the PBC-enforced phase constant
of the fast wave
. Therefore, the real part of (5) offers redundant informa-
tion, which is useful only as a confirmation that the correct phase
constant has been implemented. The imaginary part of (5) pro-
vides the attenuation constant
of the simulated
mode.
The choice of the points
and in (5) deserves special at-
tention as a consequence of the fact that the computational do-
main is the unit cell of a periodic structure and (5) is supposed to
provide an effective complex propagation constant per unit cell.
Therefore, the distance between these points should be as close
to one period
as possible. If is the cell size along the axis
of periodicity, the maximum distance between two nonboundary
Fig. 7. Field sampling points one spatial period apart from each other for the
calculation of the complex propagation constant with (5).
Fig. 8. One-dimensional leaky-wave structure treated as a phased array (under
the condition
d
).
points within the domain of one unit cell is equal to .
Evidently, this distance tends to
as the cell size tends to
zero. Alternatively, one can extend the computational domain
to include two unit cells, enforcing twice the phase difference
per unit cell at its boundaries (Fig. 7). Then, the choice of two
points exactly one period apart is straightforward, and we dis-
pense with the relevant limitation on the cell size. This approach
leads to a significant reduction in the computational cost of pe-
riodic structure modeling without compromising accuracy.
Finally, the determination of the real part of the propagation
constant as a function of frequency is pursued in the same way
as in [21].
B. Radiation Pattern Extraction
Since the modeling of NRI-based LWAs is the motivating ap-
plication of this study, the computation of radiation patterns with
the proposed FDTD analysis is also negotiated. Noting that the
unit cells of such structures are electrically small (the unit cells
of the NRI grid of [4] are less than one-tenth of a wavelength),
the phased-array approach invoked in [26] and [27] is applicable
here as well. Under this approach, each cell of the antenna is
assumed to be a point source, fed with a current of the form
(Fig. 8), where . Then, the ra-
diation pattern of the antenna is calculated using the following
expression for the corresponding normalized array factor:
(6)
KOKKINOS et al.: PERIODIC FDTD ANALYSIS OF LEAKY-WAVE STRUCTURES AND APPLICATIONS TO ANALYSIS OF NRI LWAs 1623
Fig. 9. Computational domain for the periodic FDTD analysis of the metal-
strip-loaded dielectric LWA of [28] and [29].
where is the wavelength of operation and is the angle from
the axis of propagation. Our analysis derives the
,
of an infinitely periodic structure, while (6) is evaluated for a
finite array. The number of elements
in (6) is appropriately
chosen for the computed radiation pattern to converge (typically,
is sufficient).
IV. N
UMERICAL
RESULTS:VALIDATION
A. Metal-Strip-Loaded Dielectric LWA
In this section, the proposed method is validated by compar-
ison to previously published results. A case study of interest
is the metal-strip-loaded dielectric LWA, initially proposed by
Klohn et al. [28] and then analyzed in [29]. The latter reference
follows the standard FDTD methodology of simulating a trun-
cated periodic structure composed of 29 unit cells, in order to ap-
proximate the complex propagation constant of the infinite peri-
odic one at an operating frequency of 80 GHz. On the contrary,
this study uses the compact computational domain of Fig. 9.
The dimensions involved with this geometry are as follows: the
spatial periodicity
of the structure (along the direction) is
2.5 mm. The width and the height of the dielectric substrate sup-
porting the metal strips are
mm and mm,
respectively, and the dielectric constant of the substrate is
. The width of the metal strips is . The strips are mod-
eled as perfect electric conductors. Two unit cells are included
along the
direction, while the sampling points for the electric
field are chosen to be exactly one spatial period apart. In all other
directions, the antenna is interfaced with free space, which is
terminated in Mur’s first-order absorbing boundary conditions.
The absorbing boundaries are 15 cells away from the air–dielec-
tric interfaces. In total, a mesh of
Yee’s cells is used. The computational domain is excited with
a Gabor pulse of a bandwidth of 10 GHz around the operating
frequency. This pulse excites an electric field which is polarized
parallel to the metal strips and has a spatial profile of the form
, corresponding to the
mode.
In Fig. 10, the time-domain evolution of the electric field at
a sampling point inside the dielectric substrate is shown for the
case
. This value corresponds to a leaky-wave
mode supported at 80 GHz. During the transient stage of the first
Fig. 10. Time-domain evolution of the electric field within the unit cell of the
metal-strip-loaded dielectric LWA of [28] and [29] for
=k
=
0
0
:
218
(fast-
wave mode).
Fig. 11. Complex propagation constants for the metal-strip-loaded dielectric
LWA of [28] and [29] calculated with the proposed method (periodic FDTD).
5000 time steps, the signal essentially follows the source excita-
tion. At this stage, both evanescent and propagating modes exist
in the structure. Eventually, the simulation reaches the steady
state, during which the frequency of the supported mode be-
comes dominant and a traveling wave propagates in the dielec-
tric. During this steady-state period, the amplitude of the electric
field decreases exponentially with time.
In Fig. 11, the FDTD results for the real and imaginary parts
of the complex propagation constant are presented. These re-
sults are in excellent agreement with those of [29]. Moreover, in
Fig. 12, the attenuation constant
obtained from the periodic
FDTD analysis is compared with the values of
obtained by
simulating finite versions of this LWA, with a variable number
of unit cells. For the latter simulations, Ansoft’s HFSS is used.
It is observed that
converges to the value of the infinite struc-
ture, when at least seven unit cells are employed. On the other
hand, the proposed technique relies on PBCs, thus saving cells
within the working volume and execution time.
![Fig. 11. Complex propagation constants for the metal-strip-loaded dielectric LWA of [28] and [29] calculated with the proposed method (periodic FDTD).](/figures/fig-11-complex-propagation-constants-for-the-metal-strip-1p2vctb6.webp)
![Fig. 14. Convergence of the attenuation constant (at 30 GHz) of the leaky CPW-based slot antenna array of [30], with the number of the simulated unit cells, when the antenna is analyzed as a truncated structure (with Agilent’s Momentum), compared with the proposed method (periodic FDTD).](/figures/fig-14-convergence-of-the-attenuation-constant-at-30-ghz-of-39quddxz.webp)
![Fig. 12. Convergence of the attenuation constant of the metal-strip-loaded dielectric LWA of [28] and [29], with the number of the simulated unit cells, when the antenna is analyzed as a truncated structure (with Ansoft’s HFSS), compared with the proposed method (periodic FDTD).](/figures/fig-12-convergence-of-the-attenuation-constant-of-the-metal-iq3idcjb.webp)
![Fig. 13. Computational domain for the periodic FDTD analysis of the leaky CPW-based slot antenna array of [30].](/figures/fig-13-computational-domain-for-the-periodic-fdtd-analysis-8taewh33.webp)
![Fig. 15. Top view of the 2-D loaded NRI-TL grid of [4].](/figures/fig-15-top-view-of-the-2-d-loaded-nri-tl-grid-of-4-1wxpcpgq.webp)
![Fig. 27. Zoomed-in view of the Brillouin diagram of the 2-D loaded NRI-TL grid of [4] operating as an LWA with loading elements L = 10 nH and C = 1:0 pF (closed stopband case) calculated through the periodic analysis of the equivalent circuit of Fig. 23. The zero slope of the dispersion curve at = 0 suggests that v ( = 0) = 0.](/figures/fig-27-zoomed-in-view-of-the-brillouin-diagram-of-the-2-d-e720vp6i.webp)