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 co-fired 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 ?pile-up? 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 Sol-gel 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 H-field 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 E-field 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 barium-to-strontium 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 Curie-Weiss constant fcnm Cutoff Frequency for the nth horizontal and mth vertical mode k Wave Number p(r) Charge Density Psat Saturation Polarization T0 Curie-Weiss 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. Multi-path 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 ultra-short 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 micro-fabrication 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 Sol-Gel 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, loss-less 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 ?pile-up? 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 inter-symbol 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 cross-sectional 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 non-trivial 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 ?catch-up? and ?pile-up? 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 one-dimensional 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 self-focusing and self-trapping of light beams, where energy collapses in both space and time [38,39]. If sufficient incident power is provided, the self-focusing 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 non-transverse 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 n-half 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 one-half 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 half-wavelength 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 one-half 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 non-propagating 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 quasi-infinite 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 multi-mode 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 E-field 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 quasi-TEM 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 E-probe 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 line-to-ground gap and tapered via pattern as shown in Figure 2.21. A simulated K-band 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 (24-38 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 Curie-Weiss law as given by equation (2.35), where Cc is the material-dependent Curie constant, T is the absolute temperature in kelvin and T0 is the Curie-Weiss temperature. ? = ?r?0 = ?0 + Cc|T ?T 0| (2.35) In these materials, the relative permittivity is a maximum at the Curie-Weiss temper- ature, although the Curie-Weiss 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 small-signal 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: Aband-passfilterbasedoncircularferroelectric 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 Curie-Weiss 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 Curie-Weiss 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 Curie-Weiss temperature (T0) show a small enough hysteresis loop to be considered paraelectric. Numerous studies [69?71] have confirmed this with thin-film 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 barium-to-strontium 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 low-loss 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 low-field 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 E|E| |E|?Esat ?Psat E|E| |E|?Esat (??1E+??3E3 +???) E|E| 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) E|E| (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 Curie-Weiss constant, T0 is the Curie-Weiss 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 non-ideal 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_steps-1],2,1); %input and ouput taper pop_data(4,p)=randi([3 50],1,1); %layer half-width pop_data(5:6,p)=randi([1 gene_steps-1],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,gen-1) pop_data_temp(:,1)=gen_best(:,gen-1); 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_ini-1) [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_exp-1)]); chr2=dec2bin([pop_data_temp(1:N,par2(p));(2^gene_exp-1)]); 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 half-thickness 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?,gen-1)) drawnow fprintf(?Generation %1.0i complete, reused %3.0i solutions\n?,... gen,(pop_ini-sum(gen_sol_new))); if gen>10 && gen_best(N+1,gen)==gen_best(N+1,gen-10) 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; %Non-linear 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.854e-12; %free space permittivity mu0=4*pi*1e-7; %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/(1-0.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=832e-6; %SIW hieght L_pad=0.75e-2; %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_pl-1,... (lay_start(l)+NL_reduc(i)):(lay_stop(l)-... NL_reduc(i)))=mat(l); end elseif i>=(SIW_stop-NL_taper_out(l)) NL_reduc(i)=round((lay_steps(l)-round(lay_steps(l)*... ((SIW_stop-i)/NL_taper_out(l))))/2); if NL_reduc(i)< (lay_hthick(l)) media(i,y_bound_inner_mi+1:y_bound_inner_pl-1,... (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_pl-1,... (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_mats-1)=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.x-i+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_x-1)+1; a.Ex(i)=a.max_x*(1-rho_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_x-1)+1; a.Ex(i)=a.max_x*(1-rho_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.y-j+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_y-1)+1; a.Ey(j)=a.max_y*(1-rho_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_y-1)+1; a.Ey(j)=a.max_y*(1-rho_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.z-k+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_z-1)+1; a.Ez(k)=a.max_z*(1-rho_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_z-1)+1; a.Ez(k)=a.max_z*(1-rho_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.x-i+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_x-1)+1; a.Hx(i)=H_adj*a.max_x*(1-rho_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_x-1)+1; a.Hx(i)=H_adj*a.max_x*(1-rho_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.y-1 %y negative boundary case rho_h(j)=((pml.y-j+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_y-1)+1; a.Hy(j)=H_adj*a.max_y*(1-rho_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_y-1)+1; a.Hy(j)=H_adj*a.max_y*(1-rho_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.z-1 %z negative boundary case rho_h(k)=((pml.z-k+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_z-1)+1; a.Hz(k)=H_adj*a.max_z*(1-rho_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_z-1)+1; a.Hz(k)=H_adj*a.max_z*(1-rho_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_samp-1)*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*(t-t0)).*exp(-((t-t0).^2./tau.^2)); %Gaussian modulated cosine pulse .9 cell ave fact.m %cell average factors sigEfX=0.25*(sigmaE(1:x_stop-1,2:y_stop-1,2:z_stop-1)+... sigmaE(1:x_stop-1,1:y_stop-2,2:z_stop-1)+... sigmaE(1:x_stop-1,2:y_stop-1,1:z_stop-2)+... sigmaE(1:x_stop-1,1:y_stop-2,1:z_stop-2)); sigEfY=.25*(sigmaE(2:x_stop-1,1:y_stop-1,2:z_stop-1)+... sigmaE(1:x_stop-2,1:y_stop-1,2:z_stop-1)+... sigmaE(2:x_stop-1,1:y_stop-1,1:z_stop-2)+... sigmaE(1:x_stop-2,1:y_stop-1,1:z_stop-2)); sigEfZ=.25*(sigmaE(2:x_stop-1,2:y_stop-1,1:z_stop-1)+... sigmaE(2:x_stop-1,1:y_stop-2,1:z_stop-1)+... sigmaE(1:x_stop-2,2:y_stop-1,1:z_stop-1)+... sigmaE(1:x_stop-2,1:y_stop-2,1:z_stop-1)); 146 sigMfX=(2*sigmaH(2:x_stop,1:y_stop-1,1:z_stop-1).*... sigmaH(1:x_stop-1,1:y_stop-1,1:z_stop-1))./... (sigmaH(2:x_stop,1:y_stop-1,1:z_stop-1)+... sigmaH(1:x_stop-1,1:y_stop-1,1:z_stop-1)); sigMfY=(2*sigmaH(1:x_stop-1,2:y_stop,1:z_stop-1).*... sigmaH(1:x_stop-1,1:y_stop-1,1:z_stop-1))./... (sigmaH(1:x_stop-1,2:y_stop,1:z_stop-1)+... sigmaH(1:x_stop-1,1:y_stop-1,1:z_stop-1)); sigMfZ=(2*sigmaH(1:x_stop-1,1:y_stop-1,2:z_stop).*... sigmaH(1:x_stop-1,1:y_stop-1,1:z_stop-1))./... (sigmaH(1:x_stop-1,1:y_stop-1,2:z_stop)+... sigmaH(1:x_stop-1,1:y_stop-1,1:z_stop-1)); mufX=(2*mu_eff(2:x_stop,1:y_stop-1,1:z_stop-1).*... mu_eff(1:x_stop-1,1:y_stop-1,1:z_stop-1))./... (mu_eff(2:x_stop,1:y_stop-1,1:z_stop-1)+... mu_eff(1:x_stop-1,1:y_stop-1,1:z_stop-1)); mufY=(2*mu_eff(1:x_stop-1,2:y_stop,1:z_stop-1).*... mu_eff(1:x_stop-1,1:y_stop-1,1:z_stop-1))./... (mu_eff(1:x_stop-1,2:y_stop,1:z_stop-1)+... mu_eff(1:x_stop-1,1:y_stop-1,1:z_stop-1)); mufZ=(2*mu_eff(1:x_stop-1,1:y_stop-1,2:z_stop).*... mu_eff(1:x_stop-1,1:y_stop-1,1:z_stop-1))./... 147 (mu_eff(1:x_stop-1,1:y_stop-1,2:z_stop)+... mu_eff(1:x_stop-1,1:y_stop-1,1:z_stop-1)); bexY1=repmat(be_x(2:pml.x).?,[1,y_stop-1,z_stop-2]); cexY1=repmat(ce_x(2:pml.x).?,[1,y_stop-1,z_stop-2]); bexY2=repmat(be_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-1,z_stop-2]); cexY2=repmat(ce_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-1,z_stop-2]); bexZ1=repmat(be_x(2:pml.x).?,[1,y_stop-2,z_stop-1]); cexZ1=repmat(ce_x(2:pml.x).?,[1,y_stop-2,z_stop-1]); bexZ2=repmat(be_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-2,z_stop-1]); cexZ2=repmat(ce_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-2,z_stop-1]); bhxY1=repmat(bh_x(1:pml.x).?,[1,y_stop-1,z_stop-1]); chxY1=repmat(ch_x(1:pml.x).?,[1,y_stop-1,z_stop-1]); bhxY2=repmat(bh_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-1,z_stop-1]); chxY2=repmat(ch_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-1,z_stop-1]); bhxZ=repmat(bh_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-1,z_stop-1]); chxZ=repmat(ch_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-1,z_stop-1]); bhxZ1=repmat(bh_x(1:pml.x).?,[1,y_stop-1,z_stop-1]); chxZ1=repmat(ch_x(1:pml.x).?,[1,y_stop-1,z_stop-1]); 148 bhxZ2=repmat(bh_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-1,z_stop-1]); chxZ2=repmat(ch_x(pml.x+x_steps+1:x_stop-1).?,[1,y_stop-1,z_stop-1]); .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_stop-pml.x-5; %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(n-1),w_e_out(n-1)) 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*1e-4) %if trailing edge is weak and wider than %the minimum window width x_start=max(1,trail_edge_eval-dfs+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_cent-y_bound_inner+1:y_cent+y_bound_inner-2,... z_cent-z_bound_inner+1:z_cent+z_bound_inner-1)=... (supr_start-i)^2*0.1; sigmaH(i,y_cent-y_bound_inner+1:y_cent+y_bound_inner-2,... z_cent-z_bound_inner+1:z_cent+z_bound_inner-1)=... (supr_start-i)^2*0.1/(120*pi/sqrt(fill_permE(1)))+eps; if i==x_dr sigmaE(x_dr,(y_cent-short_segs/2:y_cent),z_dr)=... sigmae_base(3); sigmaH(x_dr,(y_cent-short_segs/2:y_cent),z_dr)=... sigmah_base(3); end end cell_ave_fact; %recompute cell average factors Ex(1:(x_start-2),:,:)=0; Ey(1:(x_start-2),:,:)=0; Ez(1:(x_start-2),:,:)=0; Hx(1:(x_start-2),:,:)=0; Hy(1:(x_start-2),:,:)=0; Hz(1:(x_start-2),:,:)=0; if int_diag fprintf(?\nTrailing Boundary Advanced?) end end sc_drive=false; if (x_start-x_end-1)<=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*1e-7) && x_endpml.x+x_steps-1 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_pl-1,... z_bound_inner_mi+1:z_bound_inner_pl-1)))); plane_evalE2=max(max(E_mag(x_end-1,... y_bound_inner_mi+1:y_bound_inner_pl-1,... z_bound_inner_mi+1:z_bound_inner_pl-1))); 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_pl-1,... lay_start(l):lay_stop(l))=... E_mag_ave(SIW_start:SIW_end,... y_bound_inner_mi+1:y_bound_inner_pl-1,... lay_start(l):lay_stop(l))./E_N(diel); xiSQ(SIW_start:SIW_end,... y_bound_inner_mi+1:y_bound_inner_pl-1,... lay_start(l):lay_stop(l))=... xi(SIW_start:SIW_end,... y_bound_inner_mi+1:y_bound_inner_pl-1,... lay_start(l):lay_stop(l)).^2; ep_eff(SIW_start:SIW_end,... y_bound_inner_mi+1:y_bound_inner_pl-1,... lay_start(l):lay_stop(l))=... ep_eff(SIW_start:SIW_end,... y_bound_inner_mi+1:y_bound_inner_pl-1,... 155 lay_start(l):lay_stop(l))+... (media(SIW_start:SIW_end,... y_bound_inner_mi+1:y_bound_inner_pl-1,... 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_pl-1,... lay_start(l):lay_stop(l))+... eta(diel)^3)+xi(SIW_start:SIW_end,... y_bound_inner_mi+1:y_bound_inner_pl-1,... lay_start(l):lay_stop(l))).^(2/3)+... (sqrt(xiSQ(SIW_start:SIW_end,... y_bound_inner_mi+1:y_bound_inner_pl-1,... lay_start(l):lay_stop(l))+... eta(diel)^3)-xi(SIW_start:SIW_end,... y_bound_inner_mi+1:y_bound_inner_pl-1,... 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_end-1,2:y_stop-1,2:z_stop-1)+... ep_eff(x_start:x_end-1,1:y_stop-2,2:z_stop-1)+... ep_eff(x_start:x_end-1,2:y_stop-1,1:z_stop-2)+... 156 ep_eff(x_start:x_end-1,1:y_stop-2,1:z_stop-2)); Ca=(1-((sigEfX(x_start:x_end-1,:,:)*dt)./(2*epf)))./... (1+((sigEfX(x_start:x_end-1,:,:)*dt)./(2*epf))); Cb=(dt/epf)./(1+((sigEfX(x_start:x_end-1,:,:)*dt)./(2*epf))); Ex(x_start:x_end-1,2:y_stop-1,2:z_stop-1)=... Ca.*Exp(x_start:x_end-1,2:y_stop-1,2:z_stop-1)+... Cb.*((Hz(x_start:x_end-1,2:y_stop-1,2:z_stop-1)-... Hz(x_start:x_end-1,1:y_stop-2,2:z_stop-1))./KEyX-... (Hy(x_start:x_end-1,2:y_stop-1,2:z_stop-1)-... Hy(x_start:x_end-1,2:y_stop-1,1:z_stop-2))./KEzX); if Ni==Ni_max %update PML?s on last iteration psi_Exy(x_start:x_end-1,2:y_stop-1,2:z_stop-1)=... beyX.*psi_Exy(x_start:x_end-1,2:y_stop-1,2:z_stop-1)+... ceyX.*(Hz(x_start:x_end-1,2:y_stop-1,2:z_stop-1)-... Hz(x_start:x_end-1,1:y_stop-2,2:z_stop-1))/delta_y; Ex(x_start:x_end-1,2:y_stop-1,2:z_stop-1)=... Ex(x_start:x_end-1,2:y_stop-1,2:z_stop-1)+... Cb.*(psi_Exy(x_start:x_end-1,2:y_stop-1,2:z_stop-1)-... psi_Exz(x_start:x_end-1,2:y_stop-1,2:z_stop-1)); end %Ey epf=0.25*(ep_eff(x_start+1:x_end-1,1:y_stop-1,2:z_stop-1)+... 157 ep_eff(x_start:x_end-2,1:y_stop-1,2:z_stop-1)... +ep_eff(x_start+1:x_end-1,1:y_stop-1,1:z_stop-2)+... ep_eff(x_start:x_end-2,1:y_stop-1,1:z_stop-2)); Ca=(1-((sigEfY(x_start:x_end-2,:,:)*dt)./(2*epf)))./... (1+((sigEfY(x_start:x_end-2,:,:)*dt)./(2*epf))); Cb=(dt/epf)./(1+((sigEfY(x_start:x_end-2,:,:)*dt)./(2*epf))); Ey(x_start+1:x_end-1,1:y_stop-1,2:z_stop-1)=... Ca.*Eyp(x_start+1:x_end-1,1:y_stop-1,2:z_stop-1)+... Cb.*((Hx(x_start+1:x_end-1,1:y_stop-1,2:z_stop-1)-... Hx(x_start+1:x_end-1,1:y_stop-1,1:z_stop-2))./KEzY-... (Hz(x_start+1:x_end-1,1:y_stop-1,2:z_stop-1)-... Hz(x_start:x_end-2,1:y_stop-1,2:z_stop-1))./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_pl-1 Ey(x_dr,y_bound_inner_mi+1:y_bound_inner_pl-1,k)=... Ey(x_dr,y_bound_inner_mi+1:y_bound_inner_pl-1,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_stop-1,1:z_stop-1)=... DaX.*Hx(x_start+1:x_end,1:y_stop-1,1:z_stop-1)+... DbX.*((Ey(x_start+1:x_end,1:y_stop-1,2:z_stop)-... Ey(x_start+1:x_end,1:y_stop-1,1:z_stop-1))./KHzX-... (Ez(x_start+1:x_end,2:y_stop,1:z_stop-1)-... Ez(x_start+1:x_end,1:y_stop-1,1:z_stop-1))./KHyX); psi_Hxy(x_start+1:x_end,1:y_stop-1,1:z_stop-1)=bhyX.... *psi_Hxy(x_start+1:x_end,1:y_stop-1,1:z_stop-1)+... chyX.*(Ez(x_start+1:x_end,2:y_stop,1:z_stop-1)-... Ez(x_start+1:x_end,1:y_stop-1,1:z_stop-1))/delta_y; Hx(x_start+1:x_end,1:y_stop-1,1:z_stop-1)=... Hx(x_start+1:x_end,1:y_stop-1,1:z_stop-1)+... DbX.*(psi_Hxz(x_start+1:x_end,1:y_stop-1,1:z_stop-1)-... 162 psi_Hxy(x_start+1:x_end,1:y_stop-1,1:z_stop-1)); % Hy Hy(x_start:x_end-1,2:y_stop,1:z_stop-1)=... DaY.*Hy(x_start:x_end-1,2:y_stop,1:z_stop-1)+... DbY.*((Ez(x_start+1:x_end,2:y_stop,1:z_stop-1)-... Ez(x_start:x_end-1,2:y_stop,1:z_stop-1))./KHxY-... (Ex(x_start:x_end-1,2:y_stop,2:z_stop)-... Ex(x_start:x_end-1,2:y_stop,1:z_stop-1))./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)))>5e-6 &&... ~pulse_passed unstable=true; fprintf(?increasing energy failure detechted?) unst_type=8; end %% frame record if mod((n+frame_step-1),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)=toc-t_hold; t_hold=toc; cells=cells+(x_end-x_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_n-ave_fact 167 Pfdiv_ave(n-ave_fact+1)=... sum(Pfdiv(n-ave_fact+1:n+ave_fact-1))/(2*ave_fact-1); Pfdiv_sc_ave(n-ave_fact+1)=... sum(Pfdiv_sc(n-ave_fact+1:n+ave_fact-1))/(2*ave_fact-1); 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_end-1,:,:)*dt)./... (2*mufX(x_start:x_end-1,:,:))))./... (1+((sigMfX(x_start:x_end-1,:,:)*dt)./... (2*mufX(x_start:x_end-1,:,:)))); DbX=(dt./mufX(x_start:x_end-1,:,:))./... (1+((sigMfX(x_start:x_end-1,:,:).*dt)./... (2*mufX(x_start:x_end-1,:,:)))); KHzX=repmat(shiftdim(K.Hz(1:z_stop-1).*deltaz_H(1:z_stop-1).?,-1),... [x_end-1-(x_start-1),y_stop-1,1]); 168 KHyX=repmat(K.Hy(1:y_stop-1).*deltay(1:y_stop-1).?,[x_end-1-(x_start-1),... 1,z_stop-1]); DaY=(1-((sigMfY(x_start:x_end-1,:,:)*dt)./... (2*mufY(x_start:x_end-1,:,:))))./... (1+((sigMfY(x_start:x_end-1,:,:)*dt)./... (2*mufY(x_start:x_end-1,:,:)))); DbY=(dt./mufY(x_start:x_end-1,:,:))./... (1+((sigMfY(x_start:x_end-1,:,:).*dt)./... (2*mufY(x_start:x_end-1,:,:)))); KHxY=repmat((K.Hx(x_start:x_end-1).*deltax_H(x_start:x_end-1).?).?,... [1,y_stop-1,z_stop-1]); KHzY=repmat(shiftdim(K.Hz(1:z_stop-1).*deltaz_H(1:z_stop-1).?,-1),... [x_end-1-(x_start-1),y_stop-1,1]); DaZ=(1-((sigMfZ(x_start:x_end-1,:,:)*dt)./... (2*mufZ(x_start:x_end-1,:,:))))./... (1+((sigMfZ(x_start:x_end-1,:,:)*dt)./... (2*mufZ(x_start:x_end-1,:,:)))); DbZ=(dt./mufZ(x_start:x_end-1,:,:))./... (1+((sigMfZ(x_start:x_end-1,:,:).*dt)./... (2*mufZ(x_start:x_end-1,:,:)))); KHyZ=repmat(K.Hy(1:y_stop-1).*deltay(1:y_stop-1).?,... [x_end-1-(x_start-1),1, z_stop-1]); KHxZ=repmat((K.Hx(x_start:x_end-1).*deltax_H(x_start:x_end-1).?).?,... [1,y_stop-1,z_stop-1]); 169 beyX=repmat(be_y(1,2:y_stop-1),[x_end-x_start,1,z_stop-2]); ceyX=repmat(ce_y(1,2:y_stop-1),[x_end-x_start,1,z_stop-2]); beyZ=repmat(be_y(1,2:y_stop-1),[x_end-x_start-1,1,z_stop-1]); ceyZ=repmat(ce_y(1,2:y_stop-1),[x_end-x_start-1,1,z_stop-1]); bhyX=repmat(bh_y(1,1:y_stop-1),[x_end-x_start,1,z_stop-1]); chyX=repmat(ch_y(1,1:y_stop-1),[x_end-x_start,1,z_stop-1]); bhyZ=repmat(bh_y(1,1:y_stop-1),[x_end-x_start,1,z_stop-1]); chyZ=repmat(ch_y(1,1:y_stop-1),[x_end-x_start,1,z_stop-1]); %E field update factors KEyX=repmat(K.Ey(2:y_stop-1).*deltay(2:y_stop-1).?,[x_end-1-(x_start-1),... 1,z_stop-2]); KEzX=repmat(shiftdim(K.Ez(2:z_stop-1).*deltaz(2:z_stop-1).?,-1),... [x_end-1-(x_start-1),y_stop-2,1]); KEzY=repmat(shiftdim(K.Ez(2:z_stop-1).*deltaz(2:z_stop-1).?,-1),... [x_end-2-(x_start-1),y_stop-1,1]); KExY=repmat((K.Ex(x_start+1:x_end-1).*deltax(x_start+1:x_end-1).?).?,... [1,y_stop-1,z_stop-2]); KExZ=repmat((K.Ex(x_start+1:x_end-1).*deltax(x_start+1:x_end-1).?).?,... [1,y_stop-2,z_stop-1]); KEyZ=repmat(K.Ey(2:y_stop-1).*deltay(2:y_stop-1).?,[x_end-2-(x_start-1),... 170 1,z_stop-1]); .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_start-2:lay_stop+2),... % prop_vect(SIW_start:SIW_stop),... % squeeze(Res(1,n).EPF(SIW_start:SIW_stop,1,... % lay_start-2: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 z-y 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?,T1loc-x_dr),... sprintf(?Output: @%2.0i cells from Source?,... T2loc-x_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?,T1loc-x_dr),... sprintf(?@%2.0i cells from Source?,T2loc-x_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 X-band 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 high-pass 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