Optimized Nonlinear Substrate Integrated Waveguide for Pulse Compression
by
Byron Caudle
A dissertation submitted to the Graduate Faculty of
Auburn University
in partial fulfillment of the
requirements for the Degree of
Doctor of Philosophy
Auburn, Alabama
May 5, 2013
Keywords: FDTD, NLTL, Ferroelectrics, Pulse Compression
Copyright 2013 by Byron Caudle
Approved by
Michael Baginski, Chair, Associate Professor of Electrical and Computer Engineering
Stuart Wentworth, Associate Professor of Electrical and Computer Engineering
Lloyd Riggs, Professor of Electrical and Computer Engineering
Michael Hamilton, Assistant Professor of Electrical and Computer Engineering
Hulya Kirkici, Professor of Electrical and Computer Engineering
George Flowers, Dean of the Graduate School and Professor of Mechanical Engineering
Abstract
The time compression and associated frequency broadening of electromagnetic pulses
has numerous applications in communication and radar systems. In this work, a new type of
pulse compression device based on nonlinear ferroelectric materials in a substrate integrated
waveguide (SIW) is proposed and simulated, with low temperature cofired ceramic (LTCC)
fabrication compatibility considered. The ferroelectric material Barium Strontium Titanate
(BaxSrx?1TiO3) commonly used with tunable microwave components, is placed in vertical
slabs within the SIW and driven with an input pulse into the nonlinear polarization region.
The nonlinearity causes field amplitude dependent propagation velocity, resulting in a ten
dency for energy to ?pileup? or compress in time. The fields within the SIW are simulated
with a Finite Difference Time Domain (FDTD) method, and a Genetic Algorithm (GA) finds
the optimal material configuration that maximizes pulse compression. Pulse compression of
58% is shown by simulation to be possible with the proposed design.
ii
Acknowledgments
I would like to thank my parents, Alan and the late Ellen Caudle, for their encourage
ment and support of graduate education. I would also like to thank my wife, Kali, for her
great care and encouragement during the last two years of my Ph.D. program. Finally, I
would like to thank my advisor, Dr. Micheal Baginski and my committee for their time and
insight into the many challenging problems encountered.
iii
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
List of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1 Nonlinear Transmission Lines . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1.1 Lumped Element Designs . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1.2 Ferrite Shock Lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.1.3 Optics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Rectangular Waveguides . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.1 Waveguide Properties and Modes of Operation . . . . . . . . . . . . . 13
2.2.2 Feed Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2.3 Substrate Integrated Waveguides . . . . . . . . . . . . . . . . . . . . 22
2.3 Ferroelectric Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.1 Physical Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.2 Properties and Dependencies . . . . . . . . . . . . . . . . . . . . . . . 32
2.3.3 Loss Mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.3.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.3.5 Properties of Barium Strontium Titanate . . . . . . . . . . . . . . . . 37
2.3.6 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.3.7 Characterization and Parameter Extraction . . . . . . . . . . . . . . 45
iv
3 Principle Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.1 Theory of Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2 Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.3 Feed Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.4 Models and Initial Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4 Computational Electromagnetics Method . . . . . . . . . . . . . . . . . . . . . . 56
4.1 Finite Difference Time Domain Method and Maxwell?s Equations . . . . . . 56
4.2 FDTD Parameter Considerations . . . . . . . . . . . . . . . . . . . . . . . . 62
4.3 Absorbing Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.4 Ferroelectric Material Modeling . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.5 Source Pulse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.6 Computer Implementation and Requirements . . . . . . . . . . . . . . . . . . 72
4.7 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5 Genetic Algorithm Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.1 Process Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.2 Genes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.3 Initial Population and Parallel Processing . . . . . . . . . . . . . . . . . . . . 79
5.4 Forward Solver and Fitness Function . . . . . . . . . . . . . . . . . . . . . . 80
5.5 Ranking and Retention . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.6 Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.7 Crossover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.8 Mutation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.9 Termination Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
6 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.1 Simulation Cycle 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.2 Simulation Cycle 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
7 Fabrication Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
v
7.1 LTCC Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
7.2 Solgel Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.3 Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
8 Basic Fabricated SIW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
9 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
Appendix A: MATLAB Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
.1 GA.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
.2 FDTD main.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
.3 phys const.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
.4 E nonlin.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
.5 Control para.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
.6 UpdCo.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
.7 space ini.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
.8 sc build.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
.9 cell ave fact.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
.10 FDTD time loop.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
.11 win ups.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
.12 res plot.m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
Appendix B: FDTD Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
vi
List of Figures
2.1 Combination of dispersion and NLTL compression to create solitons . . . . . . . 4
2.2 Lumped element ladder network NLTL . . . . . . . . . . . . . . . . . . . . . . . 6
2.3 Input and output voltages for the BST filled parallel plate waveguide of Branch
and Smith [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Termination resistor used by Branch and Smith [1] for BST filled parallel pate
waveguide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5 Magnetic flux density vs. applied magnetic field, sinusoidal excitation . . . . . . 9
2.6 Coaxial ferrite shock line . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.7 Ferrite loaded shock line with axial bias [2] . . . . . . . . . . . . . . . . . . . . . 10
2.8 Top view of rectangular waveguide and reflecting wavefront vectors . . . . . . . 12
2.9 Rectangular waveguide and coordinate convention used in this work . . . . . . . 12
2.10 Electric field vectors for the TE10 mode . . . . . . . . . . . . . . . . . . . . . . 14
2.11 Top view of wave fronts at three frequencies . . . . . . . . . . . . . . . . . . . . 18
2.12 Step response of WR90 waveguide to a 40 ps rise yime step [3] . . . . . . . . . . 20
2.13 Electric field probe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.14 Magnetic loop or Hfield probe . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
vii
2.15 Standard substrate integrated waveguide structure and dimension conventions . 23
2.16 Leakage loss characteristics of substrate integrated waveguides [4] . . . . . . . . 24
2.17 Acceptable operating region for substrate integrated waveguides [4] . . . . . . . 25
2.18 Common transmission line technologies and Efield distributions . . . . . . . . . 26
2.19 Stripline with shorting posts for isolation . . . . . . . . . . . . . . . . . . . . . . 26
2.20 SIW probe feed with stripline tuning pad [5] . . . . . . . . . . . . . . . . . . . . 28
2.21 CPW to SIW transition [6] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.22 Microstrip to SIW transition [7] . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.23 Hysteresis loop formed by polarization verses applied sinusoidal electric field . . 32
2.24 Electric flux density D vs. electric field E . . . . . . . . . . . . . . . . . . . . . 33
2.25 Curie point temperature vs. barium concentration [8] . . . . . . . . . . . . . . . 38
2.26 Dielectric permittivity vs. temperature for various bariumtostrontium ratios [9] 38
2.27 Permittivity and tunability vs. temperature for BST x = 0.8 [10] . . . . . . . . 39
2.28 Hysteron model response to a sinusoidal input . . . . . . . . . . . . . . . . . . . 43
2.29 Example result of Preisach hysteresis model with 1275 hysterons . . . . . . . . . 44
2.30 Hysteresis bridge circuit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.1 Simulation Results: ax component of the Poynting vector, down the guide center 48
3.2 Rise time analysis: Minimum rise time as found 4.01 cm from the source plane . 49
viii
3.3 Cross section the nonlinear SIW . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4 SIW design with tapered end BST layers . . . . . . . . . . . . . . . . . . . . . . 52
3.5 Theoretical material set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.1 Spatial interleaving of the FDTD method . . . . . . . . . . . . . . . . . . . . . 57
4.2 Nonuniform FDTD lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.3 Cosine modulated Gaussian source pulse . . . . . . . . . . . . . . . . . . . . . . 71
4.4 Field distribution for TE10 mode source . . . . . . . . . . . . . . . . . . . . . . 71
4.5 Moving computational window components . . . . . . . . . . . . . . . . . . . . 73
5.1 Basic genetic algorithm flow chart . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.2 Desgin parameters controlled by the GA . . . . . . . . . . . . . . . . . . . . . . 78
5.3 Detailed flow chart for genetic algorithm . . . . . . . . . . . . . . . . . . . . . . 84
6.1 Simulation cycle 1: Best fitness function for each generation . . . . . . . . . . . 86
6.2 Simulation cycle 1: Optimized structure, top view . . . . . . . . . . . . . . . . . 87
6.3 Simulation cycle 1: Structure detailed top view . . . . . . . . . . . . . . . . . . 88
6.4 Simulation cycle 1: Input and output flowing power . . . . . . . . . . . . . . . . 89
6.5 Simulation cycle 1: Input and output flowing power spectrum . . . . . . . . . . 89
6.6 Simulation cycle 1: Total energy within the waveguide for all time . . . . . . . . 90
6.7 Simulation cycle 1: Planar cut of poynting vector . . . . . . . . . . . . . . . . . 92
ix
6.8 Simulation cycle 2: Optimized structure, top view . . . . . . . . . . . . . . . . . 96
6.9 Simulation cycle 2: Planar cut of poynting vector at SIW center (vertical) . . . 97
6.10 Simulation cycle 2: Input and output flowing power . . . . . . . . . . . . . . . . 98
6.11 Simulation cycle 2: Output rise time calculation detail . . . . . . . . . . . . . . 98
6.12 Simulation cycle 2: Input and output flowing power spectrum . . . . . . . . . . 99
8.1 Standard SIW and feed structure layout for 8 layers of LTCC . . . . . . . . . . 105
8.2 Simulation of fabricated standard SIW: S11 vs. frequency . . . . . . . . . . . . . 106
8.3 Simulation of fabricated standard SIW: S21 vs. frequency . . . . . . . . . . . . . 106
8.4 Fabricated standard SIW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
8.5 Fabricated standard SIW: S11 vs. frequency . . . . . . . . . . . . . . . . . . . . 107
8.6 Fabricated standard SIW: S21 vs. frequency . . . . . . . . . . . . . . . . . . . . 108
1 Code flow chart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
2 Band limited source and the associated voltage spectrum . . . . . . . . . . . . . 176
3 Wide band source and the associated voltage spectrum . . . . . . . . . . . . . . 177
4 Source electric field vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178
5 Input and output power resulting from a band limited source . . . . . . . . . . 179
6 Input and output power resulting from a wide band source . . . . . . . . . . . . 180
7 Input and output power spectrum resulting from a wide band source . . . . . . 181
x
List of Tables
2.1 Summary of prior NLTL design results . . . . . . . . . . . . . . . . . . . . . . . 5
3.1 Fixed design parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.1 FDTD simulation parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.2 Maximum values for PML parameters, based on experimental findings in [11] . . 67
5.1 Optimization parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.1 Simulation cycle 1: Genetic algorithm simulation results . . . . . . . . . . . . . 87
6.2 Simulation cycle 2: Genetic algorithm simulation results . . . . . . . . . . . . . 94
7.1 Properties of DuPontR? GreenTapeR? 9k7 [12] . . . . . . . . . . . . . . . . . . . 101
xi
List of Abbreviations
? Real Part of Propagation Constant (loss term)
?s Spatial Step
?t Time Step
? Intrinsic Impedance
? Propagation Constant
?g Guided Wavelength
B Magnetic Flux Density [Wb/m2]
D Electric Flux Density [C/m2]
E Electric Field [V/m]
H Magnetic Field [A/m]
J Current Density [A/m]
Jc Conduction Current Density [A/m]
P Polarization Density [C/m2]
? Magnetic Permeability [Webers per meter]
?0 Vacuum Magnetic Permeability = 4pi?10?7 [Webers per meter]
?r Relative Magnetic Permeability
? Conductivity [Siemens per meter]
xii
?r Rise Time
? Dielectric Permittivity [Farads per meter]
?0 Vacuum Dielectric Permittivity = 8.854?10?12 [Farads per meter]
?r Relative Dielectric Permittivity
c Speed of Light in Vacuum = 3?108 m/s
Cc CurieWeiss constant
fcnm Cutoff Frequency for the nth horizontal and mth vertical mode
k Wave Number
p(r) Charge Density
Psat Saturation Polarization
T0 CurieWeiss temperature
vg Group Velocity
vp Phase Velocity
Z0 Transmission Line Characteristic Inpedance
BST Barium Strontium Titanate (BaxStx?1TiO3)
CPML Convolutional Perfectly Matched Layer
CPW Coplanar Waveguide
FDTD Finite Difference Time Domain
GA Genetic Algorithm
LTCC Low Temperature Cofired Ceramic
xiii
NLTL Nonlinear Transmission Line
RAM Random Access Memory
SIW Substrate Integrated Waveguide
TE Transverse Electric
TEM Transverse Electric and Magnetic
TM Transverse Magnetic
xiv
Chapter 1
Introduction
The generation and transmission of short duration, broadband electromagnetic pulses
has been of significant interest to the communications and defense industries for at least
the last decade. Multipath fading, selective absorption, and multiuser interference make
narrow bandwidth (f < 5MHz) communication signals especially difficult to use in indoor
environments. Ultra wideband communications signals and radar pulses have significant
immunity to these frequency dependent effects and thus offer a distinct advantages for indoor
communications. On a larger scale, wideband radar pulses are much more difficult to jam
and largely defeat radar absorbing materials, when compared to narrow bandwidth pulses.
Pulses with ultrashort rise times are also used as trigger signals for larger circuits.
The generation, transmission and radiation of such wideband pulses has historically
been challenging due to frequency dependent dispersion and narrow bandwidth tuning cir
cuits. All transmission lines experience frequency dispersion due to parasitic reactance, thus
deviating from a true linear response to a wideband signal. These nonlinear effects can be
exploited through the creation of specialized structures known as Nonlinear Transmission
Lines (NLTLs), to change a waveform?s shape during propagation. Much work has gone into
designing nonlinear transmission lines that maximize pulse compression for the production
of solitary waves also known as solitons. A historical perspective and review of NLTL designs
is given in Chapter 2.
This works exploits the nonlinear amplitude response of the ferroelectric material Bar
ium Strontium Titanate (BaxStx?1TiO3) (BST) to cause spatial and temporal compression.
While most prior NLTL designs focused on lumped nonlinear elements in ladder networks or
1
ferrite shock lines, this work examines the use of a new type of NLTL formed by a rectan
gular waveguide. Waveguides offer an attractive alternative to microstrip traces and coaxial
cable transmission lines due to better isolation and higher power handling. The recent
development of substrate integrated waveguides (SIWs) has demonstrated that waveguide
structures can be formed at the circuit board and microfabrication level. The primary goal
of this work is the combination of classical NLTL and SIW concepts to form an optimized
pulse compression system. The design is simulated with the Finite Difference Time Domain
(FDTD) method and a Genetic Algorithm (GA) optimally selects material options and di
mensions. The design considers fabrication in Low Temperature Cofired Ceramic (LTCC)
with a SolGel process used to fill trenches (vias in the LTCC) with ferroelectric material.
A review of waveguides and ferroelectric materials necessary for development of the design
is given in Chapter 3.
2
Chapter 2
Background
2.1 Nonlinear Transmission Lines
All transmission lines have nonlinear effects that cause less than perfect, linear, lossless
propagation of electromagnetic energy. In particular, parasitic inductance and capacitance,
along with frequency dependence of the dielectric permittivity ? and magnetic permeability
?, causes waveform components to propagate with velocity as function of frequency. The
widerangeoffrequencycomponentsthatmakeupashortdurationGaussianpulsewillspread
out, or disperse, as they propagate down a real transmission line. This frequency dispersion
effect usually broadens the pulse time signature. The term ?Nonlinear Transmission Line?
refers to the intentional exploitation of nonlinear effects to achieve a specific change in a
transmitted wave, as first studied by Landauer [13]. Most NLTL designs are based on a goal
of pulse compression caused by superposition of portions of the input pulse in time. This
is accomplished by changing the propagation velocity of different portions of the pulse such
that they will ?pileup? at the same instant in time. The compressed pulse has a greater
peak power over a shorter period of time in comparison to the input pulse. As the pulse is
compressed in time, simple Fourier analysis gives an associated bandwidth broadening of the
spectral density of the pulse. The compressive effect is commonly measured as the reduction
in pulse rise time (?r) from input to output, as given by percent compression:
Compression = ?rinitial ??rfinal?
rinitial
?100; (2.1)
Some NLTL designs specifically seek to maximize the frequency spectrum broadening
effect of pulse compression. Third order and higher harmonics can be generated at high
3
Figure 2.1: Combination of dispersion and NLTL compression to create solitons
power levels, which may prove useful at frequencies and power levels beyond what can be
produced by direct synthesis [14,15].
In 1845, Russell [16] introduced the concept of a solitary water wave. These solitary
waves are now known as solitons [17]. A pure soliton propagates without loss or change in
the wave shape. If the pulse compression of an NLTL perfectly cancels the dispersion for
a given pulse shape and frequency, a soliton can propagate, as shown by Figure 2.1. Some
NLTL designs, especially optical designs, are based on this concept.
Solitons have numerous applications, but the communications applications are most
notable [18]. Device interconnects at the board level and larger are commonly limited in
data rate by frequency dispersion of square pulses. This dispersion limits the minimum rise
time and the resulting broadening of pulse eventually causes intersymbol interference if the
data rate is pushed too high. A reduced dispersion NLTL can thus be used to enable device
interconnects with higher data rates than traditional microstrip or stripline transmission
lines [19].
4
Ref NLTL Type Compression % Study Type notes
[20] Lumped Element
Nonlinear Capacitor
80% Simulated and Fab
ricated

[21] Lumped Element
Nonlinear Capacitor
82% Fabricated Nonlinear ferroelec
tric capacitors
[1] Nonlinear parallel
plate waveguide
92% Fabricated BST based nonlin
earity
[22] Ferrite shock line 98% Simulated and Fab
ricated
Uses axial magnetic
bias
Table 2.1: Summary of prior NLTL design results
Most NLTL designs can be placed in one of three categories: lumped element transmis
sion lines, ferrite shock lines, and nonlinear optical systems. Each of these categories are
discussed in the following sections. Table 2.1 summarizes the performance of several NLTL
designs.
2.1.1 Lumped Element Designs
One method of NLTL design takes advantage of either the nonlinear capacitance (C) vs.
voltage (V) curve or the nonlinear inductance (L) vs. current (I) curve of discrete specialized
capacitors or inductors respectively [23,24]. Of particular interest is the nonlinear C ?V
curve of the depletion region capacitance of a reverse bias diode [25?28]. As the depletion
region expands with increasing reverse bias (V) (limited by carrier mobility and relaxation),
the capacitance formed by the depletion interface planes will decrease as given by (2.2) and
(2.3), where A is the diode crosssectional area.
C(V) = ?0?rAd(V) (2.2)
d(V) ? 1V (2.3)
A network can be formed using these nonlinear capacitors in combination with the
inherent inductance of a microstrip transmission line, as shown in Figure 2.2. This ladder
network design has a phase velocity (vp) inversely proportional to the capacitance of each
5
Figure 2.2: Lumped element ladder network NLTL
section, as given by (2.4). Thus the frequency components of the portion of the pulse
having the largest amplitude propagate fastest (typically near the transients center). This
results in the central region of the pulse, near the peak, eventually superimposing on the
front edge, causing the front edge to appear spatially compressed with increased magnitude.
Compression islimited for anygiven designbythemaximum changeofthenonlinear elements
(?C) with reverse bias, the impedance matching of each stage to the next and the overall
impedance matching of the NLTL to the source and load impedance.
vp = 1radicalbigC(V)L(I) (2.4)
The design choices for the capacitors and inductors (or transmission line lengths) for
maximum compression is a nontrivial process; however Afshari and Hajimiri [20] demon
strated greater compression can be achieved by gradually scaling the design parameters along
the line length. A simulated and fabricated metal oxide semiconductor (MOS) varactor based
NLTL achieved 80% rise time compression.
Wilson, Turner and Smith [21] demonstrated a lumped element NLTL using nonlinear
ferroelectric filled capacitors (general ferroelectric properties are discussed in Section 2.3).
Experimental results for a 15 section network yielded 82% rise time compression. Zhao et
al. [29], recently developed simulation models that show good agreement with the fabricated
experiments using ferroelectric capacitors. Related to the prior work [21], Branch and Smith
6
Figure 2.3: Input and output voltages for the BST filled parallel plate waveguide of Branch
and Smith [1]
[1], fabricated and tested a NLTL using a parallel plate waveguide filled with BST, which
achieved at least 92% rise time compression, as shown in Figure 2.3. The challenge of creating
a matched load for the low impedance BST (? 1?2 ?) was met by terminating the parallel
plates into a copper sulfate solution, as shown by Figure 2.4. While the performance of this
NLTL design is remarkable, the implementation is impractical.
While most designs have used either nonlinear inductors or capacitors, Bostick, Nardi
and Zucker [30] have shown that hybrid lines with both nonlinear capacitors and nonlinear
inductors have greater energy storage potential. In 2012, Kuek, Liew, Schamiloglu, and
Rossi [31] examined the use of such a hybrid NLTL for the generation of microwave pulses.
When the nonlinear elements were scaled exponentially, simulations demonstrated better
matching to 50 ? resistive loads and higher amplitude output in comparison to designs
using only one type of nonlinear element.
7
Figure 2.4: Termination resistor used by Branch and Smith [1] for BST filled parallel pate
waveguide
2.1.2 Ferrite Shock Lines
Another NLTL design takes advantage of the nonlinear magnetization (B) vs. magnetic
Field (H) relationship of ferromagnetic materials. For these materials in the linear region
(small magnetic field values), the slope of the (B) vs. (H) curve represents the magnetic
permeability (?) as given by (2.5). For stronger magnetic fields, the magnetization asymp
totically approaches a near constant saturation value, proportional only to the free space
magnetic permeability ?0. The permeability for ferromagnetic materials smoothly changes
from the weak field value of ?H=0 to ?0 with increasing applied field, as shown in Figure
2.5. Given the propagation velocity of a wave through any medium as shown in (2.6) [32],
in a ferromagnetic material this velocity is nonlinear, as shown in (2.7). Thus portions of
the pulse with higher amplitude will ?catchup? and ?pileup? on the portions with smaller
amplitude, resulting in a shock wave effect [33].
8
Figure 2.5: Magnetic flux density vs. applied magnetic field, sinusoidal excitation
Figure 2.6: Coaxial ferrite shock line
B = ?H (2.5)
vp = 1??? (2.6)
vp(H) = 1radicalbig??(H) (2.7)
In 1999, Dolan [2] described and simulated a ferrite shock line based on coaxial line
with ferrite material between the center and outer conductors, as shown in Figure 2.6. Axial
magnetic bias was used to place the ferrite into the nonlinear region of the (B) vs. (H) curve
prior to the firing of the an input pulse. The test setup for such a line is shown in Figure
2.7. Sanders et al. [22] built a similar ferrite line using axial bias that reduced the 35 ns rise
time input pulse to a rise time of only 700 ps, for a rise time compression of 98%.
9
Figure 2.7: Ferrite loaded shock line with axial bias [2]
Norgard and Curry [34] presented a simulation study on an NLTL similar to ferrite shock
line using ferroelectric material instead. While this work was a very preliminary study using
onedimensional modeling, the results suggest that a coaxial line loaded with ferroelectric
materials would be successful for high power RF pulse generation.
2.1.3 Optics
Most nonlinear optics studies are concerned with the production and transmission of
solitons. Entire textbooks are devoted to the study of optical solitons [35?37]. While nonlin
ear optics and NLTL?s at microwave frequencies generally developed along separate research
lines, it is important in the context of this work to note some key similarities. Many nonlin
ear fiber optic lines utilize the Kerr nonlinearity to cancel frequency dependent (?chromatic?)
dispersion. The Kerr nonlinearity also causes selffocusing and selftrapping of light beams,
where energy collapses in both space and time [38,39]. If sufficient incident power is provided,
the selffocusing effect will continue without bound until microscopic scattering, absorption,
or breakdown results in the beam spitting into small high intensity filaments of light [40].
The intensity of this effect can even destroy the crystalline structure of some materials, as
demonstrated by the laser ablation of sapphire due to self focusing effects of air [19].
The optical Kerr effect is based on the same polarization nonlinearity as discussed in
Section 2.3 and is frequently modeled using a similar power series model. Though the effects
10
are related, the practical materials and dependencies due to temperature, frequency, etc.
are different for optical frequencies. Numerous studies [41?44] have used the FDTD method
for the study of nonlinear optical systems. The resulting detailed published material of
integrating nonlinear effects [45] made the FDTD method the clear choice for simulations to
be done in this work.
2.2 Rectangular Waveguides
Any structure that can direct electromagnetic energy in a desired direction can be
called a waveguide. In the early development of electric telegraph, simple two conductor
transmission lines were used to convey electrical current from one location to another. These
simple lines and other transmission lines types, such as coaxial cable, were very successful
in carrying communications signals. The discovery and use of electromagnetic waves for
location finding, now known as radio detection and ranging or radar, drove the need to
contain and carry high power electromagnetic energy. As the military needs leading up
to World War II demanded higher power radar systems, simple transmission lines became
impractical in terms of power loss, power handling and unintentional coupling of signals.
The solution came from using the theory and analysis of guided waves as first studied in
the late 19th century [46]. Because electromagnetic waves reflect off of a metal surface with
minimal loss, waves injected into a hollow metal tube will reflect between the walls and thus
be carried down the tube, as shown by Figure 2.8. Exploiting this concept began with the
novel use of rectangular architectural tubing to carry radio waves, but would soon become
the system of rectangular waveguides that we know today [47]. The structure of rectangular
waveguide along with the coordinate convention used for this work is detailed in Figure 2.9.
The following sections details the theory, properties and modes of rectangular waveguides.
11
Figure 2.8: Top view of rectangular waveguide and reflecting wavefront vectors
Figure 2.9: Rectangular waveguide and coordinate convention used in this work
12
2.2.1 Waveguide Properties and Modes of Operation
Basic Theory
In two conductor transmission lines, waves propagate with both electric and magnetic
field components transverse to the direction of propagation, thus known as Transverse Elec
tric and Magnetic (TEM) mode. The TEM mode satisfies Ampere?s circuit law (2.8) as long
as two conductors are present, but in a hollow waveguide, the absence of a conductor inside
the guide implies the conduction current term Jc is zero. If closed loops of magnetic field
circulate in a transverse plane (Transverse Magnetic or TM mode), as required by Gauss?s
law (?? D = 0), then (2.8) requires that there must be electric flux D and thus electric
field passing through that plane. Therefore the TM mode will have electric field terms in
the direction of propagation. Initial enforcement of a transverse electric field and applica
tion of (2.8) demonstrates a requirement for nontransverse magnetic field, giving rise to the
Transverse Electric or TE mode.
contintegraldisplay
H?dL =
integraldisplay
Jc ?dS+ ??t
integraldisplay
D?dS (2.8)
Defining the possible fields within a rectangular waveguide begins with the enforcement
of boundary conditions on Maxwell?s equations for electromagnetic waves. The waveguide
walls are initially assumed to be perfect conductors (conductivity ? = ?). Applying Fara
day?s law and Ampere?s law in integral form (enclosing the boundary between the interior of
the waveguide and the walls), gives the resulting boundary conditions that the electric field
tangent to the boundary (2.9) and magnetic field normal to the boundary (2.10) must be
zero.
Etan = 0 (2.9)
Hnorm = 0 (2.10)
13
Figure 2.10: Electric field vectors for the TE10 mode
Modes of Operation
The simplest field distribution compliant to the boundary conditions is the nhalf sinu
soidal electric field distribution across the width (?z) of the waveguide as given by the wave
number in (2.11) and the field in (2.12). The number of onehalf sinusoid variations in the
horizontal (?z) and vertical (?y) dimensions is used to identify the mode as TEnm. The first
mode, n = 1,m = 0, is referred to as lowest order mode, written as TE10 and shown in
Figure 2.10.
kc = npia (2.11)
Ey = E0 sin(kcz)e??x (2.12)
The mode distribution gives rise to the frequency  geometry dependence of waveguides,
because only electromagnetic energy with halfwavelength equal to or shorter than the waveg
uide width divided bynwill propagate in thenth mode. This maximum wavelength is known
as the cutoff wavelength as given by (2.13). The phenomenon is more commonly referred to
by the minimum frequency that will propagate, know as the cutoff frequency (2.14), where
vp is the intrinsic phase velocity within the waveguide.
14
?c10 = 2an (2.13)
fc10 = nvp2a (2.14)
Energy with frequencies higher than the cutoff frequency will also propagate; however
at double the cutoff frequency of the first mode (TE10), the next higher mode (TE20) will
propagate. This means energy at a frequency at or above the (TE20) cutoff frequency will
propagate in both (TE10) and (TE20). In practice, it is difficult to extract energy from
a waveguide in multimode operation, so waveguides are almost always limited to single
(dominant) mode (TE10) operation. The waveguide height b is limited to onehalf the width
of the guide to insure the (TE10) mode is the next higher order mode rather than the (TE11)
mode, thus maintaining the maximum operational bandwidth of the waveguide [32]. Further
reduction of the waveguide height is possible, but results in increasing loss as discussed in
Section 2.2.1.
For frequencies at or above the cutoff frequency, the wave number as given in (2.11)
can be used to find the propagation constant ?, as given by (2.15), where k is the wave
number of the dielectric inside the waveguide (generally air). The guided wavelength ?g,
given by (2.17), characteristic impedance Z0, given by (2.18), and group velocity (velocity
of flowing energy) vg, given by (2.19), can all be extracted from the propagation constant.
The intrinsic impedance ? and intrinsic phase velocity vp are given by (2.20) and (2.21)
respectively, where ? is the magnetic permeability and ? is the dielectric permittivity of the
inside of the waveguide [32,48].
? = j? = jk
radicalBigg
1?
parenleftbiggf
c
f
parenrightbigg2
(2.15)
k = 2pif??? (2.16)
15
?g = 2pi? = ?radicalbigg
1?
parenleftBig
fc
f
parenrightBig2 (2.17)
Z0 = j??? = ?radicalbigg
1?
parenleftBig
fc
f
parenrightBig2 (2.18)
vg = ???? = 1???
radicalBigg
1?
parenleftbiggf
c
f
parenrightbigg2
(2.19)
? =
radicalbigg?
? (2.20)
vp =
radicalbigg 1
?? (2.21)
The group velocity reveals two additional important properties of waveguides. First,
propagating energy near the cutoff frequency moves very slowly. This logically follows as
well from the model of waves bouncing between the waveguide walls, where the angle of
incidence for frequencies near the cutoff frequency is nearly normal to the waveguide walls.
For this reason and the subsequently discussed high conductor loss, the practical minimum
frequency for waveguides is typically 125% of the cutoff frequency. Additionaly, since group
velocity is a function of frequency, waveguides have an additional frequency dispersion effect,
which will be discussed in detail in Section 2.2.1.
Energy injected below the lowest cutoff frequency (fc10) results in a propagation constant
(2.15) that is all real, meaning such waves will be attenuated by the factore?x. Waves in this
frequency region are said to be nonpropagating or in evanescent mode. While the attenua
tion in this mode is rapid with respect to distance from the source and generally of trivial
consequence for traditional waveguide systems, it will be important to the consideration of
transient propagation in Section 2.2.1.
16
Loss Properties
Waveguide losses are due to either dielectric losses of the material filling the waveguide
(?d) or losses caused by the finite conductivity of the waveguide walls (?c). Dielectric losses
for a dry air filled guide are negligible. If the waveguide is filled with a material other than
dry air, the dielectric loss factor ??? becomes a part of the calculation of the wave number
k and leads to an wave attenuation term of e??dx on each component of the fields, as given
by (2.22), where f is the operational frequency and ? is the intrinsic impedance. For many
materials, the loss factor is stated as the loss tangent (2.23) for a specified frequency, where
? is the frequency of interest, ?? is the real part of the permittivity and ? is the conductivity
of the material filling the interior of the waveguide.
?d = 2pif?
??
2 ?
radicalBigg
1?
parenleftbiggf
c
f
parenrightbigg2
(2.22)
tan? = ??
?? +?
??? (2.23)
In thestrictestsense, assumingimperfectconductingwalls slightly changestheboundary
conditions that are used in initial waveguide analysis and therefore results in a slightly
different propagating field. The amount of change is relatively small, which allows the loss
to be calculated using perturbational methods outlined in [32]. The process can be viewed
as the summation of conductor loss around the four waveguide walls, as caused by currents
induced from magnetic fields tangent to those walls. The conductivity of the waveguide walls
must be defined by the skin effect resistance, as given by (2.25), where the permeability ?
and conductive ? are properties of the waveguide wall material. This equation also assumes
sufficiently thick waveguide walls, such that attenuation through the walls is greater than
99%.
17
Figure 2.11: Top view of wave fronts at three frequencies
?c = Rs
b?
radicalbigg
1?
parenleftBig
fc
f
parenrightBig2
bracketleftBigg
1 + 2ba
parenleftbiggf
c
f
parenrightbigg2bracketrightBigg
(2.24)
Rs =
radicalbigg
pif?
? (2.25)
The frequency dependence that appears outside of the skin effect resistance term can
be viewed by considering the reflection based waveguide model. Energy near the cutoff
frequency will reflect off the waveguide walls many more times per unit length of propagation
than energy significantly above cutoff, as depicted by Figure 2.11. Each reflection contributes
to the overall conductor loss. It is also important to note the conductor loss factor contains
an inversely proportional dependency on the waveguide height (b). If the waveguide height
is reduced from the upper limit as discussed in Section 2.2.1, conductor loss will always
increase.
Transient Propagation
The goal of this work warrants an examination of transient propagation in standard
rectangular waveguide. Though waveguide development and use was significant during
World War II, scientific investigation of the transient response did not occur until the mid
1950?s [49]. Combining an impulse with quasiinfinite frequency content or relatively content
bounded trapezoidal pulse [50] with a frequency dependent propagation constant will clearly
18
result in distortion, but the exact output from such a pulse is not readily obvious. Various
analysis [3,51?53] have sought to find more accurate and general analytical models of the
transient response. While a complete restatement of those finding will not be presented
here, the transient response effects can be summarized by combining concepts of multimode
operation, group delay, and precursor waves.
Any injected waveform with frequency content outside of the fc10 to fc20 range will
excite waves in multiple modes. This includes evanescent modes near the source, although
the impedance mismatch (Z0 for the waveguide is largely imaginary for f 0.05 (2.29)
24
Figure 2.17: Acceptable operating region for substrate integrated waveguides [4]
25
Figure 2.18: Common transmission line technologies and Efield distributions
Figure 2.19: Stripline with shorting posts for isolation
Advantages
The transmission of radio frequency and microwave signals at the circuit board level has
historically been dominated by TEM or quasiTEM transmission lines, including microstrip,
stripline and coplanar waveguide technologies. The geometry of each of these technologies
is reviewed in Figure 2.18. SIWs offer distinct advantages over each of these technologies.
Microstrip lines concentrate fields in a narrow portion of the substrate beneath the strip, but
a portion of the fields do extend into the air above the strip and outward into the substrate.
These uncontained fields can cause crosstalk between adjacent lines (both vertically and
horizontally) if strict limits on the density of lines are not maintained. Properly designed
SIWs can be placed immediately adjacent to one another, even allowing multiple SIWs
stacked on different board layers. The density limit for SIWs is essentially only a function of
frequency and the constitutive parameters (permittivity and permeability) of the material
filling the waveguide, since this dictates the waveguide width.
26
Coplanar Waveguide (CPW) is simple from a design perspective, but it suffers similar
issues to microstrip in terms of uncontained fields and possible parasitic coupling. The close
proximity of the signal trace and ground planes results in higher local electric fields than
any of the other designs, including SIWs.
Stripline designs offer better isolation than microstrip due to the upper and lower con
ducting planes. To achieve better horizontal isolation between adjacent lines requires the
use of shorting posts between the two ground planes, as shown in Figure 2.19. Such a design
differs from SIWs only by the added center conductor, making the stripline design more
complex and costly to fabricate than a SIW.
Feed Structures
As with standard waveguides, special feed structure designs are necessary to achieve
maximum power transfer in and out of the SIW. Three recently developed feed structures
are outlined below:
Yang et al. [5] designed a microstrip to Eprobe feed similar to the standard probe feed.
The thin nature of SIW?s (height?width/2) means a simple probe will be a small fraction
of a wavelength (< 0.1?) and thus relatively inefficient due to poor impedance matching.
The addition of a stripline tuning pad greatly improves this problem by reactively loading
of the probe to create a better impedance match. While the design offers DC isolation from
the input microstrip to the SIW, fabrication is complex and the performance tends to be
narrowband with return loss of ??15dB only over a 10% bandwidth.
Building on the work of Chen and Wu [56], Taringou and Dornemann [6] have recently
designed CPW to SIW transition without a backing conductor under the CPW. The design
achieved impedance matching between the CPW and the SIW by a tapered linetoground
gap and tapered via pattern as shown in Figure 2.21. A simulated Kband SIW using this
transition for input and output achieved a return loss of less than ?20dB over the entire
band. While this design is simple to fabricate, the inherent very high localized electric
27
Figure 2.20: SIW probe feed with stripline tuning pad [5]
28
Figure 2.21: CPW to SIW transition [6]
29
Figure 2.22: Microstrip to SIW transition [7]
fields of CPW mean that the total power handling capability will be limited by the CPW,
eliminating a key advantage of SIWs.
Deslandes [7] also utilized a tapered structure to create a microstrip to SIW transition,
as shown in Figure 2.22. The design builds on an empirical mode matching (for the dominant
TE10 mode) developed by Kompa [57] to find the optimum width of the the taper to SIW
junction (w). The taper length required to create the impedance match is found using the
analysis of Lu [58] and corresponds to odd multiples of ?/4. A simulated and fabricated
design using these transition on Duroid board with a 1 cm SIW achieved a return loss of
less than ?18dB over the entire TE10 bandwidth (2438 GHz).
Additional Considerations
The typical size restriction on coplanar integrated components has limited the use of
SIWs to millimeter wavelengths when these waveguides are either air filled or filled by a
low permittivity (?< 10) substrate. The use of high permittivity ferroelectric materials,
as discussed in Section 2.3, allows the use of SIWs at much lower frequency, due to much
smaller wavelengths for a given frequency. The low loss properties of Low Temperature
Cofired Ceramic (LTCC) substrates at microwave and millimeter wave frequencies make
them ideal candidates for fabricating SIWs [59,60].
30
2.3 Ferroelectric Materials
Significant research has taken place over the last 20 years into the physics of ferroelectric
materials for microwave tunable devices. While a review of the quantum mechanics analysis
commonly used is beyond the scope of this work, a brief review of ferroelectric properties
will be given, along with a discussion of material applications, the selected materials, and
modeling methods.
2.3.1 Physical Basis
The physical basis of ferroelectric material begins with the ability of certain molecules
to form an electric dipole moment. An electric dipole occurs anytime positive and negative
charges separate within a system (molecule). These dipoles can be fixed with the lattice
of a material or they can be formed by an externally applied electric field. The density of
charges making up the dipole will be denoted as p(r) . An electric dipole moment refers to
the strength of the charge separation. The density of dipole moments can be represented by
a vector field called the polarization density P . For bulk materials, it is generally acceptable
to consider the polarization density identical to the charge density p(r) distribution. This
assumption however does not consider surface states that arise at material discontinuities.
For linear materials the polarization is a simple linear function of the applied electric field,
as given by (2.30). The electric susceptibility ? is a measure of material?s ability to be
polarized.
P = ?0?E (2.30)
The polarization verses electric field response in a ferroelectric material differs in two key
ways. First, at sufficiently high electric field values the polarization begins to asymptotically
approach a saturation limit (Psat) . Secondly, after an electric field is applied and polarization
is established, a ferroelectric material will retain some polarization when the applied electric
31
Figure 2.23: Hysteresis loop formed by polarization verses applied sinusoidal electric field
field is removed. The polarization remaining when the the electric field is reduced to zero
is known as the remnant polarization. This effect is termed ?electric hysteresis? because the
polarization state is dependent on the history of both the polarization state and the presently
applied electric field. A sinusoidal cycle of electric field produces a hysteresis loop on a plot
of P vs. E, as show in Figure 2.23.
2.3.2 Properties and Dependencies
Deriving from the case of linear materials as given by (2.30), the electric flux density is
given by (2.31) where ?r = ?+1.
D = ?0E+P = ?0?rE (2.31)
In a plot of D vs. E from (2.31), as shown in Figure 2.24, the slope of the plot represents
the relative permittivity. The slope of a plot of P vs. E will only differ from the slope of D
vs. E by one. This relationship holds even when P is nonlinear in the case of a ferroelectric
32
Figure 2.24: Electric flux density D vs. electric field E
material. For such materials the general relative permittivity is given as equation (2.32) as
a function of the applied electric field.
?r (E) = 1 + ?P?E (2.32)
Thus for low values of applied electric field, the relative permittivity will be at a max
imum. As the material approaches polarization saturation, the relative permittivity will
asymptotically approach one. Many ferroelectric materials have very large values of low field
relative permittivity ?r > 100, with large variability between materials, reaching upwards of
?r = 10,000 for polarizable polymers [61].
The tunable capacitor application as discussed in Section 2.3.4 gives rise to the ?tunabil
ity? (n) or more common ?relative tunability?(nr) figures of merit for ferroelectric materials
as given by (2.33) and (2.34). E0 is a quoted value of bias field placing the material in the
nonlinear polarization region.
33
n = ?(E0)?(0) (2.33)
nr = ?(0)??(E0)?(0) (2.34)
Temperature Dependence
Many ferroelectric materials exhibit one or more phase transition temperature points.
Many materials have an electric Curie point transition above which the remnant polarization
is zero and no hysteresis loop is formed. Materials in this state are known as paraelectric or
are said to be in paraelectric phase.
The behavior of some materials can be described by the electric form of the CurieWeiss
law as given by equation (2.35), where Cc is the materialdependent Curie constant, T is the
absolute temperature in kelvin and T0 is the CurieWeiss temperature.
? = ?r?0 = ?0 + CcT ?T
0
(2.35)
In these materials, the relative permittivity is a maximum at the CurieWeiss temper
ature, although the CurieWeiss law fails due to the singularity very near T = T0. More
advanced models such as the Landau theory must be used at these temperatures [62].
Ferroelectric materials often exhibit combinations of effects, not only as functions of
temperature, but also as function of lattice stress (piezoelectric effect), fabrication sintering
temperature (often due to lattice imperfections), and dopants. Not all of these effects are
well understood. The selection of ferroelectric materials for any given application often
involves a set of complex trade offs [62]. In many cases materials with higher values of low
field permittivity and higher tunability have high dielectric loss and greater temperature
dependence of properties [63]. This work will concentrate on the use of Barium Strontium
Titanate (BST) with the known properties and considerations outlined in Section 2.3.5.
34
2.3.3 Loss Mechanisms
The general intrinsic loss mechanism within ferroelectric materials involves complex
collisions of the electromagnetic fields (and the associated contained energy) and thermal
phonons. The exact mechanisms are combinations of multiple quantum phenomena as out
lined in [62]. Extrinsic losses, such as those caused by the motion of charge inpurites and
lattice defects, can also be significant. Losses are most commonly stated as either a dielec
tric dissipation (often as a function of frequency) (2.36) or as the loss tangent at a given
frequency (2.37).
?r = ?? ?j??? (2.36)
tan? = ?+??
??
??? (2.37)
2.3.4 Applications
The ability to control relative permittivity of ferroelectric material by an externally
applied electric field gave rise to their use for tunable microwave components. In the simplest
case, a DC voltage applied across a parallel plate capacitor filled with a ferroelectric material
will produce a variable smallsignal AC capacitance as shown by (2.38), where A is the area
of the capacitor and d is the parallel plate separation distance. Since ferroelectric materials
have a large values of relative permittivity, much larger values of capacitance can be formed
with microfabrication components than is possible with only a silicon substrate.
C(V) = ?0?r(V)Ad (2.38)
Ferroelectric materials can be used in place of standard dielectric in microwave compo
nents to form tunable components from designs that would normally be narrow band. Many
devices rely on variable phase velocity through a biased ferroelectric, as given by (2.39) .
35
vp(V) = cradicalbig??
0?r(V)
(2.39)
Examplesthatutilizethiseffectinclude: Abandpassfilterbasedoncircularferroelectric
filled resonator cavities with tunable center frequency [64], a phased array lens formed by
slabs of bulk ferroelectric with biasing plates between layers [65], a phase shifter formed
by a coplanar transmission lines on ferroelectric substrates, thus having a tunable electrical
length and phase shift [66]. Three significant obstacles are faced with many microwave
device applications: high bias field requirements (E? 1MV/mtypical), variable impedance
mismatch for planer components on ferroelectric substrates (since Z0 ? 1/(??r), and lattice
interface matching problems with common substrates such as silicon [62].
While most applications utilized the bias controllable relative permittivity, one sig
nification application utilizes the hysteresis effect of ferroelectric materials. The remnant
polarization that can be held by a ferroelectric capacitor can be used as memory cell for
the fabrication of Random Access Memory (RAM) . Ferroelectric RAM or FeRAM creates a
logic ?1? by driving a ferroelectric capacitor with sufficient voltage to cause a retained spon
taneous polarization. The memory cell is read by applying a sufficient field (known as the
coercive field) to remove the remnant polarization. If the cell held a logic ?1?, a detectable
current pulse will be produced during the polarization removal process. No pulse is produced
if the cell was a logic ?0?. Unlike Dynamic RAM (DRAM), FeRAM does not require periodic
refreshing since it does not depend on maintaining a capacitor voltage which decays in time.
The primary disadvantage of FeRAM is the low memory density compared to DRAM. To
maximize the detectable current and reliability of the memory, FeRAM designs seek mate
rials with large hysteresis loops that are nearly square [67]. This is in contrast with most
microwave ferroelectric devices which work best with minimal hysteresis.
36
2.3.5 Properties of Barium Strontium Titanate
Barium Strontium TitanateBaxSrx?1TiO3 is a crystalline solid with a Perovskite struc
ture. Perovskites have a crystalline form like that calcium titanium oxide (CaTiO3), a nat
urally occurring mineral named after mineralogist L.A. Perovski. These structures take the
form ABX3, where A and B are cations of different sizes, and X is an anion, most commonly
oxygen, that bonds to both cations. The oxygen atoms are found on the face centers of the
crystals. The ferroelectric effect is common to almost all Perovskite materials [68]. While
Barium Titanate and Strontium Titanate have been used for their ferroelectric properties
independently, the ability to combine these compounds in varying fraction allows for a blend
ing of their properties. Because the Curie point of Barium Titanate and Strontium Titanate
are so vastly different (120?C vs. ?269?C respectively), blending the two allows for widely
varying control of the composite material?s Curie point, as shown in Figure 2.25 [8]. For
any given BST ratio, the permittivity approximately follows the CurieWeiss law, as shown
by Figure 2.26 [9]. Wei and Yao [19] found that both the permittivity and the tunability
are reduced when moving away from the CurieWeiss temperature (BST with x = 0.8), as
shown in Figure 2.27. For a fixed temperature, slight changes in BST ratio allows control
over material tunability.
The hysteresis properties of Barium Titanate and Strontium Titanate are also different.
Combining the two materials allows for control of the hysteresis loop. While hysteresis
effects can never be completely eliminated, BST used at an operating temperature at or
above its CurieWeiss temperature (T0) show a small enough hysteresis loop to be considered
paraelectric. Numerous studies [69?71] have confirmed this with thinfilm BST (x = 0.5)
used at room temperature. Building on the advantage for variable BST ratio, a technique has
been developed to form continuously graded BST films. Tests with tunable capacitors showed
60% tunability and much lower temperture dependance compared to uniform composition
films [72].
37
Figure 2.25: Curie point temperature vs. barium concentration [8]
Figure 2.26: Dielectric permittivity vs. temperature for various bariumtostrontium ratios
[9]
38
Figure 2.27: Permittivity and tunability vs. temperature for BST x = 0.8 [10]
While dielectric loss in BST can be attributed to the previously discussed mechanisms,
the exact loss for BST in any given application is difficult to predict, due to variability as a
function of BST ratio (x), layer thickness, sintering temperature, surface states, strain, etc.
However, some general trends can be observed. Most loss tangent measurements of BST
report values in the range 10?2 ?10?3, which is 1 to 2 orders of magnitude greater than the
loss tangent of many polymer lowloss microwave substrates (TeflonR?: tan? ? 5 ? 10?4).
High loss has been a common constraining factor from using bulk BST as tunable substrate.
Thin films (t < 5?m) of BST have a greater loss tangent than bulk forms, but for many
applications, the total loss over the film thickness is still acceptable [62]. The relatively poor
loss tangent can be improved by approximately an order of magnitude with the addition of
dopants such as magnesium oxide [73] or alumina [74]; however, the lowfield permittivity
and tunability are also reduced.
2.3.6 Modeling
Design and simulation of ferroelectric devices demands accurate and computationally
efficient models of the both polarization saturation and hysteresis effects. The following
discussion will present two methods of modeling the saturation effect and one common
method of modeling hysteresis. In all cases the models are macroscopic in nature and assume
the polarization follows the direction of the applied electric field.
39
Power Series
A piecewise power series expansion can represent the polarization as given by (2.40).
This formulation is piecewise in that a constant value of polarization saturation is used at
and above the point where an increasing electric field reaches saturation (Esat). A similar
situation occurs for the negative field sense. For low values of the electric field E, the first
term dominates (where the weak field permittivity is given as ?r = ?1 + 1). Because the
polarization for most materials is an odd function, the even order terms are zero. While many
terms can be used, in practice, only the third order term is typically used for computational
efficiency.
P =
??
????
??
???
????
+Psat EE E?Esat
?Psat EE E?Esat
(??1E+??3E3 +???) EE Otherwise
(2.40)
Hyperbolic Tangent
The double asymptotic nature of polarization saturation gives rise to considering hy
perbolic trigonometric functions exhibiting similar properties as polarization models. The
simplest model of this type is given in the (2.41) where Psat is the polarization saturation
limit and Escale is a scaling factor given by (2.42). ?r is the weak field permittivity. This
model works well for a relatively small tuning range; however, most materials approach po
larization saturation (with increasing electric field) at a slower rate than the model predicts.
P = Psat tanh(EscaleE) EE (2.41)
Escale = ?0(?r ?1)P
sat
(2.42)
40
Vendik?s Model
Vendik and Zubko [75] developed an expansion of the power series model for polarization
that can be used to directly compute the effective permittivity as a function bias field and
temperature. The model is given by (2.43)  (2.47), where E0 is the bias field, T is the
operating temperature, C is the CurieWeiss constant, T0 is the CurieWeiss temperature,
EN is the bias field scale factor and TV is a temperature adjustment factor.
?(T,E0) = ?00bracketleftBig
(?2 +?2)1/2 +?
bracketrightBig2/3
+
bracketleftBig
(?2 +?2)1/2 ??
bracketrightBig2/3
??
(2.43)
?00 = CT
0
(2.44)
? = ?00??0 (2.45)
? = E0E
n
(2.46)
? = TV?
0C
?
?
radicalBigg
1
16 +
parenleftbigg T
TV
parenrightbigg2
? T0T
V
?
? (2.47)
Vendik and Zubko [8] subsequently developed an empirical model specific for BST as a
function of the BST ratio (x), as given by (2.48)  (2.51), and the constant TV = 175[K].
T0(x) = 42 +430.37x?95.95x2 (2.48)
C(x) = parenleftbig0.86 +1.1x2parenrightbig?105 (2.49)
41
?00(x) = C(x)T
0(x)
(2.50)
EN(x) = 8.4parenleftBig
?0 (3?00(x))3/2
parenrightBig (2.51)
Preisach Hysteron Model
Ferenc (Franz) Preisach first suggested a nonideal relay or hysteron as a model for
hysteresis. This model is defined by (2.52) where k is the output state from the prior
evaluation of the input (assuming y = 0 at time zero). The response of the hysteron to a
sinusoidal input (x) is shown in Figure 2.28
y =
??
????
??
???
????
1 x??
0 x??
k ?0
tic
gens=50; %maximum number of generations
pop_ini=64; %inital population
max_l=309+1;
N=6; %number of genes
mats=10;%number of distinct material types
gene_exp=6; %number of bits per gene
gene_steps=2^gene_exp; %number of steps per gene
direct_select=floor(pop_ini*0.1);
%fraction of population directly copied
cross_prob=0.6; %Probability of crossover
mute_rate=0.02; %mutation rate
fit_list=zeros((N+1),1); %Initialize solution database
reused_sols=0; %Counter for the number of reused solution
pop_data=zeros(N+1,pop_ini); %the ?+1? index is the fitness value
%initial population fill
for p=1:pop_ini
pop_data(1,p)=randi([1 mats],1,1); %material ID
pop_data(2:3,p)=randi([1 gene_steps1],2,1); %input and ouput taper
pop_data(4,p)=randi([3 50],1,1); %layer halfwidth
pop_data(5:6,p)=randi([1 gene_steps1],2,1); %added central length
123
%and voltage step
end
%%
for gen=1:gens
%fitness function evaluation
parfor p=1:pop_ini
gen_sol_new(p)=true;
for ps=2:size(fit_list,2) %search for a prior computed solution
if all(pop_data(1:N,p)==fit_list(1:N,ps))
reused_sols=reused_sols+1;
gen_sol_new(p)=false;
Pfit(p)=fit_list(N+1,ps);
end
end
if gen_sol_new(p)
fprintf(?\nComputing solution %1.0i in gen %1.0i\n?,...
p,gen)
mat_vect=pop_data(1,p)+5; %material index vector
t_in=pop_data(2,p); %Input taper
t_out=pop_data(3,p); %Output taper
lay_half_thick=pop_data(4,p);
lay_length=pop_data(5,p);
V_scale=pop_data(6,p);
[Pfit(p),unstable(p)]=...
FDTD_main(mat_vect,t_in,t_out,lay_half_thick,...
lay_length,V_scale);
124
end
end
for p=1:pop_ini %add generated fitness solutions
%to fit_list database
if gen_sol_new(p)
fit_list=[fit_list [(pop_data(1:N,p)); Pfit(p)]];
end
end
pop_data(N+1,:)=Pfit;
[fit_order,pop_order]=sort(pop_data(N+1,:),2,?descend?);
%sort fitness function
for p=1:pop_ini %sort individuals by fitness
pop_data_temp(:,p)=pop_data(:,pop_order(p));
end
%Elitism Operation
if gen>1 && pop_data_temp(N+1,1)<=gen_best(N+1,gen1)
pop_data_temp(:,1)=gen_best(:,gen1);
end
if any(unstable) %Report unstable solutions
unst_count=sum(unstable);
fprintf(?FDTD instability failure on %i NEW solutions?,...
125
unst_count)
end
%Roulette wheel parameters
wheel_size=sum((pop_data_temp(N+1,:)));
wheel_space=cumsum((pop_data_temp(N+1,:)));
wheel_prob=wheel_space./wheel_size;
%direct copy to next generation
pop_data(:,(1:direct_select))=pop_data_temp(:,(1:direct_select));
for p=direct_select+1:2:(pop_ini1)
[p1temp,par1(p)]=min(abs(rand(1)wheel_prob));
[p2temp,par2(p)]=min(abs(rand(1)wheel_prob));
chr1=dec2bin([pop_data_temp(1:N,par1(p));(2^gene_exp1)]);
chr2=dec2bin([pop_data_temp(1:N,par2(p));(2^gene_exp1)]);
if rand(1)1 %prevent mutation of layer id bits
if rand(1)<=mute_rate
if child_chr(g,b)==?0?
child_chr(g,b)=?1?;
else
child_chr(g,b)=?0?;
end
end
end
end
end
pop_data(1:N,p)=bin2dec(child_chr(1:N,:));
pop_data(1,p)=min((mats)*ones(1,1),(pop_data(1,p)));
%clamp to maximum material id #
pop_data(4,p)=min(50,(pop_data(4,p)));
%clamp to maximum layer halfthickness
127
pop_data(4,p)=max(3,(pop_data(4,p)));
%clamp layer thickness to minimum value
pop_data(1:N,p)=max(ones(N,1),(pop_data(1:N,p)));
%clamp all to minimum value of 1
end
bar(pop_data_temp(N+1,:))
title(sprintf(?Sorted fitness for generation: %3.0i?,gen1))
drawnow
fprintf(?Generation %1.0i complete, reused %3.0i solutions\n?,...
gen,(pop_inisum(gen_sol_new)));
if gen>10 && gen_best(N+1,gen)==gen_best(N+1,gen10)
fprintf(?\nSimulation halted on appearent convergence \n?)
gen_best(:,gen)
break
end
end
matlabpool close
t_end=toc/60
cl_data=clock;
data_fl=[?GA14_? num2str(cl_data(2)) ?_? num2str(cl_data(3)),...
?_? num2str(cl_data(4)) ?h? num2str(cl_data(5)) ?m?];
save(data_fl);
128
else
fprintf(?\nMATLABpool failed to start?)
end
.2 FDTD main.m
function [sol_fit,unstable]=FDTD_main(mat,NL_taper_in,NL_taper_out,...
lay_hthick,lay_length,V_scale)
int_plot=false; %controls the generation of internal test plots
int_diag=false; %controls the generation of internal diagnostics
if int_diag
tic;
end
%% Setup parameters and coefficients
phsy_const; %physical constants
max_N=5000; frame_step=30; %time steps and frame step
E_nonlin; %Nonlinear permitivitty
Control_Para; %Control Parameters
UpdCo; %Update Coefficients and CPML
space_ini; %intialize the field space
sc_build; %Source setup
cell_ave_fact; %Averaging of constituative parameters
if int_diag
t_setup=toc; %Measure setup time
end
%% Main FDTD run
sim_run=true; unstable=false; %Booleans to control simulation stop
129
%on instability or high computational error
sol_fit=0; %initialize fitness function
FDTD_time_loop;%Main FDTD time marching loop
if int_diag
t_end=toc; %run time for the simulation run
tic
res_plot
end
end
.3 phys const.m
%Physical Constants
ep0=8.854e12; %free space permittivity
mu0=4*pi*1e7; %free space permeability
c=1/sqrt(mu0*ep0); %speed of light
eta0=sqrt(mu0/ep0); %free space intrinsic impedance
.4 E nonlin.m
%Nonlinear material modeling parameters
dl_mats=10; %number of material options
BSTx=linspace(0.5,0.25,dl_mats); %setup material set BST ratio?s
T=300; %operating temperature in K
T_V=175;
for diel=1:dl_mats;
T_0(diel)=42+439.37*BSTx(diel)95.95*BSTx(diel)^2;
C(diel)=(0.86+1.1*BSTx(diel)^2)/(10^5);
ep00(diel)=C(diel)/T_0(diel);
130
alpha(diel)=(T_V/(ep0*C(diel)))*(sqrt((1/16)+(T/T_V)^2)...
(T_0(diel)/T_V));
eta(diel)=ep00(diel)*alpha(diel)*ep0;
E_N(diel)=(1/(10.7))*8.4/(ep0*(3*ep00(diel))^(3/2));
ep_0bias(diel)=1/(alpha(diel)*ep0)./10; %divided by 10 due to dopant
%reduced of permittivity
end
.5 Control para.m
%Control Parameters for 8 layer 9k7 LTCC SIW
pulse_amp=V_scale*100;
ep_min_dt=7.1; %mimimum permittivity used for time step calculation
fill_permH=1; %permeability scale factor
an_rate=0.03;%frame delay in results animation
%%parameters / dimensions / spacing
a_dim=0.008581; %SIW width
b_dim=832e6; %SIW hieght
L_pad=0.75e2; %Length of SIW on each end of the BST layer
wall_thick_z=2; %wall thickness # of segments
wall_thick_y=2; %wall thickness # of segments
f_max=40e9; %maximum frequency of interest
lam_minLTCC=(3e8/sqrt(7.1)/f_max); %minimum wavelength in LTCC
N_y_LTCC=10; %number of cells at minimum wavelength in LTCC
delta_min_LTCC=lam_minLTCC/N_y_LTCC; %largest spatial step size in LTCC
N_y=10; %minimum wavelength in BST
lam_min=(3e8/sqrt(173)/f_max); %number of cells at minimum wavelength in BST
131
delta_min=lam_min/N_y; %spatial step size in BST
delta_y=delta_min*2; %spatial step size in y
short_segs=floor(b_dim/delta_y); %number of segment for the waveguide height
deltay=ones(short_segs,1)*delta_y; %vector represent step size in y
expand_fact=1.1; %rate of expansion of the nonuniform lattice
lay_steps=2*lay_hthick; %total thickness of the BST layer
lay_start=0+wall_thick_z; lay_stop=lay_steps; %layer start and stop points
deltaz=ones(lay_steps,1)*delta_min; %vector represent step size in z
wg_width=sum(deltaz); %variable used to accumulate total structure width
delta_min_holdx=expand_fact*delta_min; %step size within loop
while wg_width1
media(i,y_bound_inner_mi+1:y_bound_inner_pl1,...
(lay_start(l)+NL_reduc(i)):(lay_stop(l)...
NL_reduc(i)))=mat(l);
end
elseif i>=(SIW_stopNL_taper_out(l))
NL_reduc(i)=round((lay_steps(l)round(lay_steps(l)*...
((SIW_stopi)/NL_taper_out(l))))/2);
if NL_reduc(i)< (lay_hthick(l))
media(i,y_bound_inner_mi+1:y_bound_inner_pl1,...
(lay_start(l)+NL_reduc(i)):(lay_stop(l)...
NL_reduc(i)))=mat(l);
end
else %full width region
NL_reduc(i)=0;
media(i,y_bound_inner_mi+1:y_bound_inner_pl1,...
(lay_start(l)+NL_reduc(i)):(lay_stop(l)...
NL_reduc(i)))=mat(l);
end
end
%material definitions
137
ep_base(1:4)=ep0;
ep_base(5)=7.1*ep0;
ep_base(6:6+dl_mats1)=ep_0bias*ep0;
mu_base(1:size(ep_base,2))=mu0;
sigmae_base=zeros(1,size(ep_base,2));
sigmae_base(3)=3.5e7;
omega_loss=10e9*2*pi;
loss_tan=0.02;
sigmae_base(6:Nmedia)=loss_tan*(ep_0bias*ep0*omega_loss);
sigmah_base=zeros(1,size(ep_base,2))+eps;
.6 UpdCo.m
%% Material parameter preallocation
sigma.Ex=zeros(1,x_stop)+eps;
sigma.Hx=zeros(1,x_stop)+eps;
sigma.Ey=zeros(1,y_stop)+eps;
sigma.Hy=zeros(1,y_stop)+eps;
sigma.Ez=zeros(1,z_stop)+eps;
sigma.Hz=zeros(1,z_stop)+eps;
a.Ex=zeros(1,x_stop)+eps;
a.Hx=zeros(1,x_stop)+eps;
a.Ey=zeros(1,y_stop)+eps;
a.Hy=zeros(1,y_stop)+eps;
a.Ez=zeros(1,z_stop)+eps;
a.Hz=zeros(1,z_stop)+eps;
138
%Averaged effective parameters
sigmaE=zeros(x_stop,y_stop,z_stop);
sigmaH=zeros(x_stop,y_stop,z_stop);
ep_eff=zeros(x_stop,y_stop,z_stop);
mu_eff=zeros(x_stop,y_stop,z_stop);
%% E sigma and K factors
for i=1:1:x_stop %Loop for x dependent values
if i <= pml.x %x negative boundary case
rho_e(i)=((pml.xi+0.75)/(pml.x));
poly_fact=rho_e(i)^m_pml;
sigma.Ex(i)=poly_fact*sigma.max_x;
K.Ex(i)=poly_fact*(K.max_x1)+1;
a.Ex(i)=a.max_x*(1rho_e(i));
elseif i > pml.x+x_steps %x positive boundary case
edge=x_steps+pml.x; %last free space cell
rho_e(i)=(i(edge+1)+0.25)/(pml.x);
poly_fact=rho_e(i)^m_pml;
sigma.Ex(i)=poly_fact*sigma.max_x;
K.Ex(i)=poly_fact*(K.max_x1)+1;
a.Ex(i)=a.max_x*(1rho_e(i));
else
K.Ex(i)=1; %central computation space
end
be_x(i)=exp((sigma.Ex(i)/K.Ex(i)+a.Ex(i))*(dt/ep0));
ce_x(i)=(sigma.Ex(i)*(be_x(i)1))/(K.Ex(i)*...
139
(sigma.Ex(i)+a.Ex(i)*K.Ex(i)));
end
clear rho_e
for j=1:1:y_stop %Loop for y dependent values
if j <= pml.y %y negative boundary case
rho_e(j)=((pml.yj+0.75)/(pml.y));
poly_fact=rho_e(j)^m_pml;
sigma.Ey(j)=poly_fact*sigma.max_y;
K.Ey(j)=poly_fact*(K.max_y1)+1;
a.Ey(j)=a.max_y*(1rho_e(j));
elseif j > pml.y+y_steps %y positive boundary case
edge=y_steps+pml.y; %first pml cell
rho_e(j)=(j(edge+1)+0.25)/(pml.y);
poly_fact=rho_e(j)^m_pml;
sigma.Ey(j)=poly_fact*sigma.max_y;
K.Ey(j)=poly_fact*(K.max_y1)+1;
a.Ey(j)=a.max_y*(1rho_e(j));
else
K.Ey(j)=1; %central computation space
end
be_y(j)=exp((sigma.Ey(j)/K.Ey(j)+a.Ey(j))*(dt/ep0));
ce_y(j)=(sigma.Ey(j)*(be_y(j)1))/(K.Ey(j)*...
(sigma.Ey(j)+a.Ey(j)*K.Ey(j)));
end
140
clear rho_e
for k=1:1:z_stop %Loop for z dependent values
if k <= pml.z %z negative boundary case
rho_e(k)=((pml.zk+0.75)/(pml.z));
poly_fact=rho_e(k)^m_pml;
sigma.Ez(k)=poly_fact*sigma.max_z;
K.Ez(k)=poly_fact*(K.max_z1)+1;
a.Ez(k)=a.max_z*(1rho_e(k));
elseif k > pml.z+z_steps %z positive boundary case
edge=pml.z+z_steps; %first pml cell
rho_e(k)=(k(edge+1)+0.25)/(pml.z);
poly_fact=rho_e(k)^m_pml;
sigma.Ez(k)=poly_fact*sigma.max_z;
K.Ez(k)=poly_fact*(K.max_z1)+1;
a.Ez(k)=a.max_z*(1rho_e(k));
else
K.Ez(k)=1; %central computation space
end
be_z(k)=exp((sigma.Ez(k)/K.Ez(k)+a.Ez(k))*(dt/ep0));
ce_z(k)=(sigma.Ez(k)*(be_z(k)1))/(K.Ez(k)*...
(sigma.Ez(k)+a.Ez(k)*K.Ez(k)));
end
%% H update half offset sigma and K factors
141
offset=0.5;
H_adj=1;
for i=1:1:x_stop %Loop for x dependent values
if i <= pml.x %x negative boundary case
rho_h(i)=((pml.xi+0.25)/(pml.x));
poly_fact=rho_h(i)^m_pml;
sigma.Hx(i)=H_adj*poly_fact*sigma.max_x;
K.Hx(i)=poly_fact*(K.max_x1)+1;
a.Hx(i)=H_adj*a.max_x*(1rho_h(i));
elseif i > pml.x+x_steps %x positive boundary case
edge=x_steps+pml.x; %last free space cell
rho_h(i)=(i(edge+1)+0.75)/(pml.x);
poly_fact=rho_h(i)^m_pml;
sigma.Hx(i)=H_adj*poly_fact*sigma.max_x;
K.Hx(i)=poly_fact*(K.max_x1)+1;
a.Hx(i)=H_adj*a.max_x*(1rho_h(i));
else
K.Hx(i)=1;
end
bh_x(i)=exp((sigma.Hx(i)/K.Hx(i)+a.Hx(i))*(dt/ep0));
ch_x(i)=(sigma.Hx(i)*(bh_x(i)1))/(K.Hx(i)*...
(sigma.Hx(i)+a.Hx(i)*K.Hx(i)));
end
clear rho_h
for j=1:1:y_stop %Loop for y dependent values
142
if j <= pml.y1 %y negative boundary case
rho_h(j)=((pml.yj+0.25)/(pml.y));
poly_fact=rho_h(j)^m_pml;
sigma.Hy(j)=H_adj*poly_fact*sigma.max_y;
K.Hy(j)=poly_fact*(K.max_y1)+1;
a.Hy(j)=H_adj*a.max_y*(1rho_h(j));
elseif j > pml.y+y_steps %y positive boundary case
edge=y_steps+pml.y; %first pml cell
rho_h(j)=(j(edge+1)+0.75)/(pml.y);
poly_fact=rho_h(j)^m_pml;
sigma.Hy(j)=H_adj*poly_fact*sigma.max_y;
K.Hy(j)=poly_fact*(K.max_y1)+1;
a.Hy(j)=H_adj*a.max_y*(1rho_h(j));
else
K.Hy(j)=1; %central computation space
end
bh_y(j)=exp((sigma.Hy(j)/K.Hy(j)+a.Hy(j))*(dt/ep0));
ch_y(j)=(sigma.Hy(j)*(bh_y(j)1))/(K.Hy(j)*...
(sigma.Hy(j)+a.Hy(j)*K.Hy(j)));
end
clear rho_h
for k=1:1:z_stop %Loop for z dependent values
if k <= pml.z1 %z negative boundary case
rho_h(k)=((pml.zk+0.25)/(pml.z));
poly_fact=rho_h(k)^m_pml;
143
sigma.Hz(k)=H_adj*poly_fact*sigma.max_z;
K.Hz(k)=poly_fact*(K.max_z1)+1;
a.Hz(k)=H_adj*a.max_z*(1rho_h(k));
elseif k > pml.z+z_steps %z positive boundary case
edge=pml.z+z_steps; %first pml cell
rho_h(k)=(k(edge+1)+0.75)/(pml.z);
poly_fact=rho_h(k)^m_pml;
sigma.Hz(k)=H_adj*poly_fact*sigma.max_z;
K.Hz(k)=poly_fact*(K.max_z1)+1;
a.Hz(k)=H_adj*a.max_z*(1rho_h(k));
else
K.Hz(k)=1;%central computation space
end
bh_z(k)=exp((sigma.Hz(k)/K.Hz(k)+a.Hz(k))*(dt/ep0));
ch_z(k)=(sigma.Hz(k)*(bh_z(k)1))/(K.Hz(k)*...
(sigma.Hz(k)+a.Hz(k)*K.Hz(k)));
end
%% Initialize some parameter arrays
for m=1:Nmedia
m_id=(media==m);
sigmaE=sigmaE+m_id*sigmae_base(m);
sigmaH=sigmaH+m_id*sigmah_base(m);
ep_eff=ep_eff+m_id*ep_base(m);
mu_eff=mu_eff+m_id*mu_base(m);
end
144
.7 space ini.m
%Intialize the field space
Ex=zeros(x_stop,y_stop,z_stop);
Hx=zeros(x_stop,y_stop,z_stop);
Ey=zeros(x_stop,y_stop,z_stop);
Hy=zeros(x_stop,y_stop,z_stop);
Ez=zeros(x_stop,y_stop,z_stop);
Hz=zeros(x_stop,y_stop,z_stop);
psi_Exy=zeros(x_stop,y_stop,z_stop);
psi_Exz=zeros(x_stop,y_stop,z_stop);
psi_Eyx=zeros(x_stop,y_stop,z_stop);
psi_Eyz=zeros(x_stop,y_stop,z_stop);
psi_Ezx=zeros(x_stop,y_stop,z_stop);
psi_Ezy=zeros(x_stop,y_stop,z_stop);
psi_Hxy=zeros(x_stop,y_stop,z_stop);
psi_Hxz=zeros(x_stop,y_stop,z_stop);
psi_Hyx=zeros(x_stop,y_stop,z_stop);
psi_Hyz=zeros(x_stop,y_stop,z_stop);
psi_Hzx=zeros(x_stop,y_stop,z_stop);
psi_Hzy=zeros(x_stop,y_stop,z_stop);
.8 sc build.m
%Source Builder
f_cut=1/(2*a_dim*sqrt(ep0*mu0*7.1*fill_permH)); %TE_10 cutoff frequency
145
f_cent=f_cut*1.5; %waveguide center band frequency
T=dt; %sample time
N_samp=50000; %maximum number of time samples
t=(0:N_samp1)*T; %time vector
tau=550*T; %Gaussian variance
t0=900*T; %time step offset
stop_n=2000; %last value of ?n? with any significant source value
Vs=pulse_amp.*cos(2*pi*f_cent*(tt0)).*exp(((tt0).^2./tau.^2));
%Gaussian modulated cosine pulse
.9 cell ave fact.m
%cell average factors
sigEfX=0.25*(sigmaE(1:x_stop1,2:y_stop1,2:z_stop1)+...
sigmaE(1:x_stop1,1:y_stop2,2:z_stop1)+...
sigmaE(1:x_stop1,2:y_stop1,1:z_stop2)+...
sigmaE(1:x_stop1,1:y_stop2,1:z_stop2));
sigEfY=.25*(sigmaE(2:x_stop1,1:y_stop1,2:z_stop1)+...
sigmaE(1:x_stop2,1:y_stop1,2:z_stop1)+...
sigmaE(2:x_stop1,1:y_stop1,1:z_stop2)+...
sigmaE(1:x_stop2,1:y_stop1,1:z_stop2));
sigEfZ=.25*(sigmaE(2:x_stop1,2:y_stop1,1:z_stop1)+...
sigmaE(2:x_stop1,1:y_stop2,1:z_stop1)+...
sigmaE(1:x_stop2,2:y_stop1,1:z_stop1)+...
sigmaE(1:x_stop2,1:y_stop2,1:z_stop1));
146
sigMfX=(2*sigmaH(2:x_stop,1:y_stop1,1:z_stop1).*...
sigmaH(1:x_stop1,1:y_stop1,1:z_stop1))./...
(sigmaH(2:x_stop,1:y_stop1,1:z_stop1)+...
sigmaH(1:x_stop1,1:y_stop1,1:z_stop1));
sigMfY=(2*sigmaH(1:x_stop1,2:y_stop,1:z_stop1).*...
sigmaH(1:x_stop1,1:y_stop1,1:z_stop1))./...
(sigmaH(1:x_stop1,2:y_stop,1:z_stop1)+...
sigmaH(1:x_stop1,1:y_stop1,1:z_stop1));
sigMfZ=(2*sigmaH(1:x_stop1,1:y_stop1,2:z_stop).*...
sigmaH(1:x_stop1,1:y_stop1,1:z_stop1))./...
(sigmaH(1:x_stop1,1:y_stop1,2:z_stop)+...
sigmaH(1:x_stop1,1:y_stop1,1:z_stop1));
mufX=(2*mu_eff(2:x_stop,1:y_stop1,1:z_stop1).*...
mu_eff(1:x_stop1,1:y_stop1,1:z_stop1))./...
(mu_eff(2:x_stop,1:y_stop1,1:z_stop1)+...
mu_eff(1:x_stop1,1:y_stop1,1:z_stop1));
mufY=(2*mu_eff(1:x_stop1,2:y_stop,1:z_stop1).*...
mu_eff(1:x_stop1,1:y_stop1,1:z_stop1))./...
(mu_eff(1:x_stop1,2:y_stop,1:z_stop1)+...
mu_eff(1:x_stop1,1:y_stop1,1:z_stop1));
mufZ=(2*mu_eff(1:x_stop1,1:y_stop1,2:z_stop).*...
mu_eff(1:x_stop1,1:y_stop1,1:z_stop1))./...
147
(mu_eff(1:x_stop1,1:y_stop1,2:z_stop)+...
mu_eff(1:x_stop1,1:y_stop1,1:z_stop1));
bexY1=repmat(be_x(2:pml.x).?,[1,y_stop1,z_stop2]);
cexY1=repmat(ce_x(2:pml.x).?,[1,y_stop1,z_stop2]);
bexY2=repmat(be_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop1,z_stop2]);
cexY2=repmat(ce_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop1,z_stop2]);
bexZ1=repmat(be_x(2:pml.x).?,[1,y_stop2,z_stop1]);
cexZ1=repmat(ce_x(2:pml.x).?,[1,y_stop2,z_stop1]);
bexZ2=repmat(be_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop2,z_stop1]);
cexZ2=repmat(ce_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop2,z_stop1]);
bhxY1=repmat(bh_x(1:pml.x).?,[1,y_stop1,z_stop1]);
chxY1=repmat(ch_x(1:pml.x).?,[1,y_stop1,z_stop1]);
bhxY2=repmat(bh_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop1,z_stop1]);
chxY2=repmat(ch_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop1,z_stop1]);
bhxZ=repmat(bh_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop1,z_stop1]);
chxZ=repmat(ch_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop1,z_stop1]);
bhxZ1=repmat(bh_x(1:pml.x).?,[1,y_stop1,z_stop1]);
chxZ1=repmat(ch_x(1:pml.x).?,[1,y_stop1,z_stop1]);
148
bhxZ2=repmat(bh_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop1,z_stop1]);
chxZ2=repmat(ch_x(pml.x+x_steps+1:x_stop1).?,[1,y_stop1,z_stop1]);
.10 FDTD time loop.m
%FDTD main time loop
%little ?n? for time step variable
frame_count=0; %initialize frame count for result recording
cells=0; %Total number of cells computed (for statistics only)
supress_region=8; %must be less the value used below x_dr+____.
plane_evalE1=0; plane_evalE2=0; x_start=1;
dfs=5; %distance from source for evaluation
trail_edge_eval=x_dr+dfs;
supr_start=trail_edge_eval; %start point of trailing edge suppression zone
if int_diag
t_hold=t_setup;
end
unstable=false; unst_type=0; %Instability stop controls
pp_ct=0; %count of time steps passed since ?pulse_passed? idicator tripped
ft_look=x_stoppml.x5; %output evaluation plane (cell #)
w_e_in(1)=0; w_e_out(1)=0; %initializae input and output energy storage
win_ups;
% P_flow=zeros(x_stop,N);
pulse_passed=false; %boolean for pulse passing through output
final_n=N; %last time step computed (updated elsewhere)
or_ind=0; %Overwrite index for storage of E_mag
%%
149
for n=1:1:N
%% Time step based controls
win_move=false; %window movement indicator
if int_diag
clc
fprintf(?%i%% complete  step %i  %i to %i  run #%i\n?,...
round(n/N*100),n,x_start,x_end,run_count)
if n>1
fprintf(?%e joules in, %e joules out?,w_e_in(n1),w_e_out(n1))
end
end
if n<(stop_n) %during time source is active
sc_drive=true;
elseif n==stop_n %at the switch off of the source
win_min=x_end; %Window is prevented from shrink below it width
%when the source switched off
elseif n>stop_n %begin the possiblity of moving the trail window edge.
if plane_evalE1<(pulse_amp/delta_y*1e4)
%if trailing edge is weak and wider than
%the minimum window width
x_start=max(1,trail_edge_evaldfs+1); %advance the trailing edge
trail_edge_eval=trail_edge_eval+1;
supr_start=supr_start+1; %advance the edge
%of the suppression region
win_move=true;
150
for i=x_start:supr_start %setup suppresion region
sigmaE(i,y_centy_bound_inner+1:y_cent+y_bound_inner2,...
z_centz_bound_inner+1:z_cent+z_bound_inner1)=...
(supr_starti)^2*0.1;
sigmaH(i,y_centy_bound_inner+1:y_cent+y_bound_inner2,...
z_centz_bound_inner+1:z_cent+z_bound_inner1)=...
(supr_starti)^2*0.1/(120*pi/sqrt(fill_permE(1)))+eps;
if i==x_dr
sigmaE(x_dr,(y_centshort_segs/2:y_cent),z_dr)=...
sigmae_base(3);
sigmaH(x_dr,(y_centshort_segs/2:y_cent),z_dr)=...
sigmah_base(3);
end
end
cell_ave_fact; %recompute cell average factors
Ex(1:(x_start2),:,:)=0;
Ey(1:(x_start2),:,:)=0;
Ez(1:(x_start2),:,:)=0;
Hx(1:(x_start2),:,:)=0;
Hy(1:(x_start2),:,:)=0;
Hz(1:(x_start2),:,:)=0;
if int_diag
fprintf(?\nTrailing Boundary Advanced?)
end
end
sc_drive=false;
if (x_startx_end1)<=win_min && int_diag;
151
%detect minimum window operation
fprintf(?\nAt or below minimum window Operation?)
end
if trail_edge_eval+10==x_end
fprintf(strcat(?\nSimulation halted on wave passing out ?,...
?of computation region\n?))
sim_run=false;
break
end
%
end
%leading window edge movement
if x_end==2*pml.x+x_steps %do nothing if window has reached max point
if int_diag
fprintf(?\nSimulation reached PML boundary?)
end
else %conditional expansion ofthe leading edge of the window
if plane_evalE2>(pulse_amp*1e7) && x_endpml.x+x_steps1
x_end=x_stop; %snap in PML as wave approaches
else
x_end=x_end+1; %otherwise gradually expand window
end
end
152
end
if plane_evalE2>50e6  isnan(plane_evalE1) %instability evaluation
unstable=true;
unst_type=2;
end
if plane_evalE1>pulse_amp*50000  isnan(plane_evalE1)
%instability evaluation
unstable=true;
unst_type=3;
end
if win_move %if windows has moveed or expanded, recacluated
%coefficient arrays
win_ups;
end
%% Field magnitude and prior value holding and instability evaluations
E_mag(1:x_end,:,:)=sqrt(Ex(1:x_end,:,:).^2+Ey(1:x_end,:,:).^2+...
Ez(1:x_end,:,:).^2);
plane_evalE1=max(max(max(E_mag(trail_edge_eval:trail_edge_eval+10,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
z_bound_inner_mi+1:z_bound_inner_pl1))));
plane_evalE2=max(max(E_mag(x_end1,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
z_bound_inner_mi+1:z_bound_inner_pl1)));
if max(max(max(E_mag)))>120e6; %breakdown voltage check
153
unstable=true;
unst_type=7;
end
E_mag_sf=20; %time averaging for E_mag used
%to calculate nonlinear permittivity
if n>E_mag_sf
or_ind=or_ind+1;
if or_ind>E_mag_sf
or_ind=1;
end
E_mag_store(SIW_start:SIW_stop,:,lay_start:lay_stop,or_ind)=...
sqrt(Ex(SIW_start:SIW_stop,:,lay_start:lay_stop).^2+...
Ey(SIW_start:SIW_stop,:,lay_start:lay_stop).^2+...
Ez(SIW_start:SIW_stop,:,lay_start:lay_stop).^2);
E_mag_ave=sum(E_mag_store,4)/E_mag_sf;
else
or_ind=or_ind+1;
E_mag_store(SIW_start:SIW_stop,:,lay_start:lay_stop,or_ind)=...
sqrt(Ex(SIW_start:SIW_stop,:,lay_start:lay_stop).^2+...
Ey(SIW_start:SIW_stop,:,lay_start:lay_stop).^2+...
Ez(SIW_start:SIW_stop,:,lay_start:lay_stop).^2);
E_mag_ave(SIW_start:SIW_stop,:,lay_start:lay_stop)=...
E_mag_store(SIW_start:SIW_stop,:,lay_start:lay_stop,n);
end
%% Update ep_eff components and E fields
154
Exp=Ex; Eyp=Ey; Ezp=Ez; %prior value holding
SIW_end=min(SIW_stop,x_end);
Ni_max=2; %number of iteration in Newton iteration
%for nonlinear materials
for Ni=1:1:Ni_max
ep_eff(SIW_start:SIW_end,:,:)=(media(SIW_start:SIW_end,:,:)<6).*...
ep_eff(SIW_start:SIW_end,:,:);
for l=1:lay_ct
if mat(l)>=6
diel=mat(l)5;
xi(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l))=...
E_mag_ave(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l))./E_N(diel);
xiSQ(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l))=...
xi(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l)).^2;
ep_eff(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l))=...
ep_eff(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
155
lay_start(l):lay_stop(l))+...
(media(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l))==...
mat(l)).*max(ep_min_dt,ep00(diel)./...
((sqrt(xiSQ(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l))+...
eta(diel)^3)+xi(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l))).^(2/3)+...
(sqrt(xiSQ(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l))+...
eta(diel)^3)xi(SIW_start:SIW_end,...
y_bound_inner_mi+1:y_bound_inner_pl1,...
lay_start(l):lay_stop(l))).^(2/3)eta(diel))).*ep0./10;
%divided by 10 due to dopant
end
end
%Ex
epf=0.25*(ep_eff(x_start:x_end1,2:y_stop1,2:z_stop1)+...
ep_eff(x_start:x_end1,1:y_stop2,2:z_stop1)+...
ep_eff(x_start:x_end1,2:y_stop1,1:z_stop2)+...
156
ep_eff(x_start:x_end1,1:y_stop2,1:z_stop2));
Ca=(1((sigEfX(x_start:x_end1,:,:)*dt)./(2*epf)))./...
(1+((sigEfX(x_start:x_end1,:,:)*dt)./(2*epf)));
Cb=(dt/epf)./(1+((sigEfX(x_start:x_end1,:,:)*dt)./(2*epf)));
Ex(x_start:x_end1,2:y_stop1,2:z_stop1)=...
Ca.*Exp(x_start:x_end1,2:y_stop1,2:z_stop1)+...
Cb.*((Hz(x_start:x_end1,2:y_stop1,2:z_stop1)...
Hz(x_start:x_end1,1:y_stop2,2:z_stop1))./KEyX...
(Hy(x_start:x_end1,2:y_stop1,2:z_stop1)...
Hy(x_start:x_end1,2:y_stop1,1:z_stop2))./KEzX);
if Ni==Ni_max %update PML?s on last iteration
psi_Exy(x_start:x_end1,2:y_stop1,2:z_stop1)=...
beyX.*psi_Exy(x_start:x_end1,2:y_stop1,2:z_stop1)+...
ceyX.*(Hz(x_start:x_end1,2:y_stop1,2:z_stop1)...
Hz(x_start:x_end1,1:y_stop2,2:z_stop1))/delta_y;
Ex(x_start:x_end1,2:y_stop1,2:z_stop1)=...
Ex(x_start:x_end1,2:y_stop1,2:z_stop1)+...
Cb.*(psi_Exy(x_start:x_end1,2:y_stop1,2:z_stop1)...
psi_Exz(x_start:x_end1,2:y_stop1,2:z_stop1));
end
%Ey
epf=0.25*(ep_eff(x_start+1:x_end1,1:y_stop1,2:z_stop1)+...
157
ep_eff(x_start:x_end2,1:y_stop1,2:z_stop1)...
+ep_eff(x_start+1:x_end1,1:y_stop1,1:z_stop2)+...
ep_eff(x_start:x_end2,1:y_stop1,1:z_stop2));
Ca=(1((sigEfY(x_start:x_end2,:,:)*dt)./(2*epf)))./...
(1+((sigEfY(x_start:x_end2,:,:)*dt)./(2*epf)));
Cb=(dt/epf)./(1+((sigEfY(x_start:x_end2,:,:)*dt)./(2*epf)));
Ey(x_start+1:x_end1,1:y_stop1,2:z_stop1)=...
Ca.*Eyp(x_start+1:x_end1,1:y_stop1,2:z_stop1)+...
Cb.*((Hx(x_start+1:x_end1,1:y_stop1,2:z_stop1)...
Hx(x_start+1:x_end1,1:y_stop1,1:z_stop2))./KEzY...
(Hz(x_start+1:x_end1,1:y_stop1,2:z_stop1)...
Hz(x_start:x_end2,1:y_stop1,2:z_stop1))./KExY);
if Ni==Ni_max %update PML?s on last iteration
if x_startE_mag_sf %recompute E_mag_ave with new values
E_mag_store(SIW_start:SIW_stop,:,...
lay_start:lay_stop,E_mag_sf)=...
sqrt(Ex(SIW_start:SIW_stop,:,lay_start:lay_stop).^2+...
Ey(SIW_start:SIW_stop,:,lay_start:lay_stop).^2+...
Ez(SIW_start:SIW_stop,:,lay_start:lay_stop).^2);
E_mag_ave=sum(E_mag_store,4)/E_mag_sf;
end
end %end nonlinear iteration loop
161
%Source Drive
if n<=stop_n
for k=z_bound_inner_mi+1:z_bound_inner_pl1
Ey(x_dr,y_bound_inner_mi+1:y_bound_inner_pl1,k)=...
Ey(x_dr,y_bound_inner_mi+1:y_bound_inner_pl1,k)+...
Vs(n)/((short_segs)*delta_y)*cos((pi/2)*sc_taper_fact(k));
end
end
%% Update H components
% Hx
Hx(x_start+1:x_end,1:y_stop1,1:z_stop1)=...
DaX.*Hx(x_start+1:x_end,1:y_stop1,1:z_stop1)+...
DbX.*((Ey(x_start+1:x_end,1:y_stop1,2:z_stop)...
Ey(x_start+1:x_end,1:y_stop1,1:z_stop1))./KHzX...
(Ez(x_start+1:x_end,2:y_stop,1:z_stop1)...
Ez(x_start+1:x_end,1:y_stop1,1:z_stop1))./KHyX);
psi_Hxy(x_start+1:x_end,1:y_stop1,1:z_stop1)=bhyX....
*psi_Hxy(x_start+1:x_end,1:y_stop1,1:z_stop1)+...
chyX.*(Ez(x_start+1:x_end,2:y_stop,1:z_stop1)...
Ez(x_start+1:x_end,1:y_stop1,1:z_stop1))/delta_y;
Hx(x_start+1:x_end,1:y_stop1,1:z_stop1)=...
Hx(x_start+1:x_end,1:y_stop1,1:z_stop1)+...
DbX.*(psi_Hxz(x_start+1:x_end,1:y_stop1,1:z_stop1)...
162
psi_Hxy(x_start+1:x_end,1:y_stop1,1:z_stop1));
% Hy
Hy(x_start:x_end1,2:y_stop,1:z_stop1)=...
DaY.*Hy(x_start:x_end1,2:y_stop,1:z_stop1)+...
DbY.*((Ez(x_start+1:x_end,2:y_stop,1:z_stop1)...
Ez(x_start:x_end1,2:y_stop,1:z_stop1))./KHxY...
(Ex(x_start:x_end1,2:y_stop,2:z_stop)...
Ex(x_start:x_end1,2:y_stop,1:z_stop1))./KHzY);
if x_start(t0/T/2) && max(Pfdiv)>max(Pfdiv_sc) && pulse_passed==0
pulse_passed=true;
165
end
if n>(t0/T/8) && max(diff(w_e(SIW_start+1:end,n)))>5e6 &&...
~pulse_passed
unstable=true;
fprintf(?increasing energy failure detechted?)
unst_type=8;
end
%% frame record
if mod((n+frame_step1),frame_step)==0 && int_diag
frame_count=frame_count+1;
Res(frame_count).Ey=Ey(:,y_store,:);
Res(frame_count).Ez=Ez(:,y_store,:);
Res(frame_count).EPF=ep_eff(:,y_store,:);
Res(frame_count).Px=Px(:,y_store,:);
plot_bounds(frame_count,:)=[x_start x_end];
end
%% time step stats
if int_diag
t_step(n)=toct_hold;
t_hold=toc;
cells=cells+(x_endx_start)*y_stop*z_stop;
end
%% Unstable halt
if any(isnan(P_flow(:,n)))
unstable=true;
unst_type=6;
end
166
if unstable
fprintf(?\n\nSimulation numerical instability, type %i\n\n?,...
unst_type)
sim_run=false;
final_n=n;
Pfit=0; Pfdiv_max=0;
break
end
if pulse_passed
pp_ct=pp_ct+1;
if pp_ct>400 %count the time steps after the pulse has past
sim_run=false;
final_n=n;
break
end
end
%% Simulation completion
if n==N %simulation has run to completion
sim_run=false;
final_n=n;
end
end
%% Solution Evaluation
if unstable==0
ave_fact=3;
for n=ave_fact:final_nave_fact
167
Pfdiv_ave(nave_fact+1)=...
sum(Pfdiv(nave_fact+1:n+ave_fact1))/(2*ave_fact1);
Pfdiv_sc_ave(nave_fact+1)=...
sum(Pfdiv_sc(nave_fact+1:n+ave_fact1))/(2*ave_fact1);
end
Pfdiv_sc_max=max(Pfdiv_sc_ave);
Pfdiv_max=max(Pfdiv_ave);
sol_fit=Pfdiv_max/Pfdiv_sc_max
end
if int_diag
fprintf(?Final simulation step ran from x cell %f to %f\n?,...
x_start,x_end);
end
.11 win ups.m
%Window dependant update factors
%H field update factors
DaX=(1((sigMfX(x_start:x_end1,:,:)*dt)./...
(2*mufX(x_start:x_end1,:,:))))./...
(1+((sigMfX(x_start:x_end1,:,:)*dt)./...
(2*mufX(x_start:x_end1,:,:))));
DbX=(dt./mufX(x_start:x_end1,:,:))./...
(1+((sigMfX(x_start:x_end1,:,:).*dt)./...
(2*mufX(x_start:x_end1,:,:))));
KHzX=repmat(shiftdim(K.Hz(1:z_stop1).*deltaz_H(1:z_stop1).?,1),...
[x_end1(x_start1),y_stop1,1]);
168
KHyX=repmat(K.Hy(1:y_stop1).*deltay(1:y_stop1).?,[x_end1(x_start1),...
1,z_stop1]);
DaY=(1((sigMfY(x_start:x_end1,:,:)*dt)./...
(2*mufY(x_start:x_end1,:,:))))./...
(1+((sigMfY(x_start:x_end1,:,:)*dt)./...
(2*mufY(x_start:x_end1,:,:))));
DbY=(dt./mufY(x_start:x_end1,:,:))./...
(1+((sigMfY(x_start:x_end1,:,:).*dt)./...
(2*mufY(x_start:x_end1,:,:))));
KHxY=repmat((K.Hx(x_start:x_end1).*deltax_H(x_start:x_end1).?).?,...
[1,y_stop1,z_stop1]);
KHzY=repmat(shiftdim(K.Hz(1:z_stop1).*deltaz_H(1:z_stop1).?,1),...
[x_end1(x_start1),y_stop1,1]);
DaZ=(1((sigMfZ(x_start:x_end1,:,:)*dt)./...
(2*mufZ(x_start:x_end1,:,:))))./...
(1+((sigMfZ(x_start:x_end1,:,:)*dt)./...
(2*mufZ(x_start:x_end1,:,:))));
DbZ=(dt./mufZ(x_start:x_end1,:,:))./...
(1+((sigMfZ(x_start:x_end1,:,:).*dt)./...
(2*mufZ(x_start:x_end1,:,:))));
KHyZ=repmat(K.Hy(1:y_stop1).*deltay(1:y_stop1).?,...
[x_end1(x_start1),1, z_stop1]);
KHxZ=repmat((K.Hx(x_start:x_end1).*deltax_H(x_start:x_end1).?).?,...
[1,y_stop1,z_stop1]);
169
beyX=repmat(be_y(1,2:y_stop1),[x_endx_start,1,z_stop2]);
ceyX=repmat(ce_y(1,2:y_stop1),[x_endx_start,1,z_stop2]);
beyZ=repmat(be_y(1,2:y_stop1),[x_endx_start1,1,z_stop1]);
ceyZ=repmat(ce_y(1,2:y_stop1),[x_endx_start1,1,z_stop1]);
bhyX=repmat(bh_y(1,1:y_stop1),[x_endx_start,1,z_stop1]);
chyX=repmat(ch_y(1,1:y_stop1),[x_endx_start,1,z_stop1]);
bhyZ=repmat(bh_y(1,1:y_stop1),[x_endx_start,1,z_stop1]);
chyZ=repmat(ch_y(1,1:y_stop1),[x_endx_start,1,z_stop1]);
%E field update factors
KEyX=repmat(K.Ey(2:y_stop1).*deltay(2:y_stop1).?,[x_end1(x_start1),...
1,z_stop2]);
KEzX=repmat(shiftdim(K.Ez(2:z_stop1).*deltaz(2:z_stop1).?,1),...
[x_end1(x_start1),y_stop2,1]);
KEzY=repmat(shiftdim(K.Ez(2:z_stop1).*deltaz(2:z_stop1).?,1),...
[x_end2(x_start1),y_stop1,1]);
KExY=repmat((K.Ex(x_start+1:x_end1).*deltax(x_start+1:x_end1).?).?,...
[1,y_stop1,z_stop2]);
KExZ=repmat((K.Ex(x_start+1:x_end1).*deltax(x_start+1:x_end1).?).?,...
[1,y_stop2,z_stop1]);
KEyZ=repmat(K.Ey(2:y_stop1).*deltay(2:y_stop1).?,[x_end2(x_start1),...
170
1,z_stop1]);
.12 res plot.m
%results plotter
FR=frame_count; %total number of frames
view(114,14)
%% Main plot
for n=1:FR
[az,el]=view;
surf(trans_vect,prop_vect(1:x_stop),(squeeze(Res(1,n).Px(:,1,:)./...
(100^2)))); scale_lim=50e7; zlim([1,1].*scale_lim);
% surf(trans_vect,prop_vect(1:x_stop),(squeeze(Res(1,n).Ey(:,1,:)...
% 0*media(:,y_cent,:)))); scale_lim=50e6; zlim([1,1].*scale_lim);
% surf(trans_vect(lay_start2:lay_stop+2),...
% prop_vect(SIW_start:SIW_stop),...
% squeeze(Res(1,n).EPF(SIW_start:SIW_stop,1,...
% lay_start2:lay_stop+2)./ep0));
% scale_lim=175; zlim([0,1].*scale_lim);
zlim([1,1].*scale_lim);
set(gca,?CLim?,[.4,1].*scale_lim);
ylabel(?Propagation (m)?)
xlabel(?Transverse (m)?)
zlabel(?E_y amplitude at the center of the guide (V/m)?)
% zlabel(?P_x (W/cm^2) @ y center ?)
shading flat
171
% colormap(gray)
grid on
title([?Time step #?,num2str(n*frame_step)])
view(az,el) %recycle current positon
drawnow
pause(an_rate)
end
%% y center cut
clf; set(gcf,?visible?,?on?);
for n=1:10:final_n
plot(prop_vect(1:size(P_flow,1)),P_flow(1:end,n));
scale_lim=6e7;
ylim([1,1].*scale_lim);
ylabel(?Total P_{flow} [W] through zy planes?)
title([?Time step #?,num2str(n)])
xlabel(?Propagation Dimension [m]?)
grid on
pause(an_rate)
end
%% P_flow at loc vs. time
clf;
T2loc=ft_look;
T1loc=x_dr+10;
172
PfT1=P_flow(T1loc,:);
PfT2=P_flow(T2loc,:);
subplot(2,1,1)
hold on
plot([1:1:size(P_flow,2)]*dt,[PfT1].?)
plot([1:1:size(PfT2,2)]*dt,[PfT2].?,?r?)
hold off
grid minor
xlabel(?Time [s]?)
ylabel(?P_{flow}?)
title(?P_{flow} vs. Time?)
legend(sprintf(?Input @%2.0i cells from Source?,T1locx_dr),...
sprintf(?Output: @%2.0i cells from Source?,...
T2locx_dr),?location?,?best?);
fs=dt;
NFFT=2^(nextpow2(size(P_flow,2))+5);
Pff1=fft(PfT1,NFFT)/size(P_flow,2);
Pff2=fft(PfT2,NFFT)/size(PfT2,2);
f_scale=Fs*linspace(0,1,NFFT);
subplot(2,1,2)
hold on
plot(f_scale(1:end)./1e9,abs([Pff1(1:end)].?))
plot(f_scale(1:end)./1e9,abs([Pff2(1:end)].?),?r?)
hold off
xlim([0 50])
xlabel(?Frequency [GHz]?)
173
ylabel(?P_{flow} (f)?)
legend(sprintf(?@%2.0i cells from Source?,T1locx_dr),...
sprintf(?@%2.0i cells from Source?,T2locx_dr));
grid minor
%% Total energy vs space
for n=1:20:size(P_flow,2)
plot([sum(P_flow(:,1:n),2)*dt, sum(P_flow,2)*dt]);
xlabel(?Spatial Step?)
ylabel(?Total Energy?)
title(sprintf(?Through time step %1.0i?,n))
drawnow
pause(an_rate)
end
174
Appendix B: FDTD Validation
The core FDTD algorithm (without nonlinearity) was validated by simulation of the
fabricated Xband structure in Chapter 8. A nonuniform grid was used with the minimum
lattice space set to ?min = ?min/60, where ?min = c/(??rLTCCfmax). The maximum fre
quency of interest was chosen as twice the center operating frequency: fmax = 20 GHz.
Simulation used either a band limited (fc10 to fc20) source (Figure 2) or a wide band Gaus
sian pulse source (Figure 3). The source was applied to the microstrip as a soft source setting
of vertical field components between the microstrip trace and the ground plane, as given in
(1), and shown in Figure 4, where b is the substrate thickness (SIW height).
Eny(xdr,y,z) = Eny(xdr,y,z) + Vs(n)b (1)
175
0 1 2 3 4 5 6
x 10?10
?1000
?500
0
500
1000
Time step
V source
(V)
0 2 4 6 8 10 12 14 16 180
5
10
15
frequency in GHz
V source
(f)
Single mode band
Figure 2: Band limited source and the associated voltage spectrum
176
0 0.2 0.4 0.6 0.8 1 1.2 1.4
x 10?10
0
200
400
600
800
1000
Time step
V source
(V)
0 2 4 6 8 10 12 14 16 180
1
2
3
frequency in GHz
V source
(f)
Single mode band
Figure 3: Wide band source and the associated voltage spectrum
177
Figure 4: Source electric field vectors
178
0 0.5 1 1.5
x 10?9
?0.2
?0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
Time [s]
P flow
normalized to maximum input
Power through input plane
Power through output plane
Output peak
at 87% of input peak
Small reflection due to transition
Figure 5: Input and output power resulting from a band limited source
When the band limited source is used, a very small reflection occurs due to the tapered
microstrip transitions, as can be seen in Figure 5. This is expected since the transition length
is a function of wavelength. The large conductivity of the assumed silver metalization results
in a small loss from input to output of only ? 0.6 dB.
179
0 0.5 1 1.5
x 10?9
?0.2
?0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
Time [s]
P flow
normalized to maximum input
Input @10 cells from Source
Output: @426 cells from Source
Large reflection
at transition
Only the small portion
of energy above
fc10 passes
Figure 6: Input and output power resulting from a wide band source
Using the wide band Gaussian pulse, significant reflections can be observed returning
through the input plane, as represented by negative values in Figure 6. The spectrum
of the input power, as shown in Figure 7, is a result of the complex combination of the
multiple reflections. A close examination of the output power spectrum reveals a deep null
at f ? 2fc10, which corresponds to the expected highpass filtering effect of the waveguide.
This explains why only a small option of the power passes from input to output.
180
0 10 20 30 40 50?0.2
0
0.2
0.4
0.6
0.8
1
1.2
Frequency [GHz]
P
flow
(f) normalized to maximum input
Input @10 cells from Source
Output: @426 cells from Source
?Rough? spectrum
due to multiple
reflections
Deep spectral null
at ~2fc10
Figure 7: Input and output power spectrum resulting from a wide band source
181