IEEE MTT-V044-I06 (1996-06)

IEEE T R A N S A C T I 0 N S

ON

MICROWAVE THEORY AND TECHNIQUES A PUBLICATION OF THE IEEE MICROWAVE THEORY AND TECHNIQUES SOCIETY

JUNE 1996

VOLUME 44

NUMBER 6

IETMAB

(ISSN 0018-9480)

[email protected]

PAPERS Stepped transformers on TEM-transmission lines - V.P. Meschanov ; I.A. Rasukova ; V.D. Tupikin Electromagnetic wave effects on microwave transistors using a full-wave time-domain model –

793

- 798

M.A. Alsunaidi ; S.M.S. Imtiaz ; S.M. El-Ghazaly 799 Finite-element method with edge elements for waveguides loaded with ferrite magnetized in arbitrary direction - Lezhu Zhou ; L.E. Davis 809 Influence of higher order modes on the measurements of complex permittivity and permeability of materials using a microstrip discontinuity P. Queffelec ; P. Gelin 816 Design and measurements of a novel subharmonically pumped millimeter-wave mixer using two single planar Schottky-barrier diodes S.D. Vogel 825 Finite-difference time-domain algorithm for solving Maxwell's equations in rotationally symmetric geometries Yinchao Chen ; R. Mittra ; P. Harms 832 Characterization of buried microstrip lines for constructing high-density microwave integrated circuits - T. Ishikawa ; E. Yamashita 840 Study of whispering gallery modes in double disk sapphire resonators - Hong Peng 848 Modeling of general constitutive relationships in SCN TLM - L.R.A.X. de Menezes ; W.J.R. Hoefer 854 High-frequency reciprocity based circuit model for the incidence of electromagnetic waves on general circuits in layered media F. Olyslager 862 Higher-order vector finite elements for tetrahedral cells - J.S. Savage ; A.F. Peterson 874 Riccati matrix differential equation formulation for the analysis of nonuniform multiple coupled microstrip lines - Jen-Tsai Kuo 880 Field simulation of dipole antennas for interstitial microwave hyperthermia - G. Schaller ; J. Erb ; R. Engelbrecht 887 CAD models for multilayered substrate interdigital capacitors - S.S. Gevorgian ; T. Martinsson ; P.L.J. Linner ; E.L. Kollberg 896 Computer-aided design and optimization of NRD-guide mode suppressors - Jifu Huang ; Ke Wu ; F. Kuroki ; T. Yoneyama 905 A new miniature magnetic field probe for measuring three-dimensional fields in planar high-frequency circuits - Yingjie Gao ; I. Wolff 911 Parameter extraction and correction for transmission lines and discontinuities using the finite-difference time-domain method M.A. Schamberger ; S. Kosanovich ; R. Mittra 919 Open-ended metallized ceramic coaxial probe for high-temperature dielectric properties measurements - S. Bringhurst ; M.F. Iskander 926 DC instability of the series connection of tunneling diodes - O. Boric-Lubecke ; Dee-Son Pan ; T. Itoh 936 Mode coupling in superconducting parallel plate resonator in a cavity with outer conductive enclosure –

- 808 - 815

Feng Gao ; M.V. Klein ; J. Kruse ; Milton Feng Coupling parameters for a side-coupled ring resonator and a microstrip line - Shih-Lin Lu ; A.M. Ferendeci The use of active traveling-wave structures in GaAs MMIC's - S.G. Ingram ; J.C. Clifton Finite-difference time-domain analysis of flip-chip interconnects with staggered bumps - H.H.M. Ghouz ; E.-B. El-Sharawy

- 952 - 956 - 960 - 963

( Continued on back cover)

944 953 956 960

- 824 - 831 - 839 - 847 - 853 - 861 - 873 - 879 - 886 - 895 - 904 - 910 - 918 - 925 - 935 - 943

Thermal management for high-power active amplifier arrays - N.J. Kolias ; R.C. Compton 963 LSE- and LSM-mode sheet impedances of thin conductors - S. Amari ; J. Bornemann 967 Accurate analysis of losses in waveguide structures by compact two-dimensional FDTD method combined with autoregressive signal analysis M. Fujii ; S. Kobayashi 970 A simple formula for the concentration of charge on a three-dimensional corner of a conductor - Yimin Zhang ; A.H. Zemanian 975 Computation of equivalent circuits of CPW discontinuities using quasi-static spectral domain method - D. Mirshekar-Syahkal 979 High-power HTS planar filters with novel back-side coupling - Zhi-Yuan Shen ; C. Wilker ; P. Pang ; C. Carter 984 An integrated SIS mixer and HEMT IF amplifier –

- 966 - 970

S. Padin ; D.P. Woody ; J.A. Stern ; H.G. LeDuc ; R. Blundell ; C.-Y.E. Tong ; M.W. Pospieszalski

- 990

(end)

987

- 975 - 979 - 984 - 986

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

791

Stepped Transformers on TEM-Transmission Lines Valery Petrovich Meschanov, Member, ZEEE, Irina Anatolyevna Rasukova, and Vladimir Dmitrievich Tupikin

Abstract-The paper presents comparative analysis of the properties of impedance stepped transformers both with monotonous and nonmonotonous step-to-step impedance variation. A miniature stepped transformer of a new structure based on a cascade of an even number of uniform transmission line sections has been synthesized. Section lengths are considerably shorter than a quarter of the central wavelength, and the section impedances alternate. The proposed transformers are the simplest to implement among the available analogs. As an example, the results of the solution of the Chebyshev approximation problem for the four-and six-section transformers of different specifications are given.

I. INTRODUCTION

T

RANSFORMERS are traditionally divided into two groups, the first comprising the transmission-line devices with a continuously tapered impedance function (nonuniform transformers), the second including those on the transmission lines with piecewise constant variation of the impedance function (stepped transformers). The latter are considerably shorter than the tapered transformers and find broad application. Stepped transformers are further divided into monotonic (with a monotonic step-to-step impedance variation) and nonmonotonic (with nonmonotonic step-to-step impedance variation). The monotonic transformers, which are to be regarded as classic, have been proposed and investigated in sufficient detail by the American scientists [ 11-[3], while the nonmonotonic ones have been described mainly by the Russian authors [41-[71. The present paper gives the comparative analysis of the properties for both monotonic and nonmonotonic stepped transformers. The results of investigation of the properties for miniature nonmonotonic stepped transformers having a new structure and characterized both by minimum length, optimum frequency characteristics, and distinguished from the available analogs by their simplicity of design are also presented.

lines to be matched, x, ( i = 1 , 2 , . . . , n ) are impedances of quarterwave sections). According to the classifications introduced in [4], the transformers under consideration refer to stepped transformers of Class I. They are antimetry devices; according to Riblet [3], the antimetry condition may be written in the form Z,Z,+1-,

=

xz

i = 1.2, . . . , n.

(1)

The main drawback of Class I transformers is their considerable length L = nX,/4, where ri is the number of transformer sections. Stepped transformers of Class I1 synthesized using m cascaded uniform transmission line sections of various lengths with alternating impedances ( m is always an even number) are shorter by a factor of 1.5-2 [5]. In a particular case the section impedances are equal to the impedances of the transmission lines to be matched [ 5 ] , [ti], and [8] [see Fig. l(b)]. The elementary two-section devices of such a type have been proposed previously [SI; but they did not find wide application on account of their narrow bandwidth. More complex multistepped units have not been investigated. Recently the Chebyshev approximation problem for the prescribed amplitude frequency characteristic for multistep superwide band impedance transformers of Class I1 has been solved by the authors [5]. The problem has been formulated as follows: to define the component values for the vector A = ( A I ,A2, . . . , A m ) allowing one to achieve

where 6' = 27rX0/X is the generalized electric variable;

81,& correspond to the lower and the upper matching band boundaries; X is the transmission line wavelength; II'(6', A)I is the modulus of the input reflection coefficient, Ir(6', A)I =

(l-l/~5!'11~z)1~2,where T11 is the element of the wave transfer matrix for the transformer. The vector A components are the normalized section lengths L, = l,/& ( i .= 1 , 2 , . . . , m ) , where 1, denotes the geometrical lengths of the sections. 11. GENERALPROPERTIES OF THE STEPPED TRANSFORMERS The solution of the corresponding approximation problems [5] has led the authors to the statement that the optimum Properties of the stepped transformers on cascade connec- Chebyshev characteristics can be provided only by the struction of n uniform transmission line sections of equal lengths ture, for which the relations (3) are true 1, = X0/4 (A, is the wave length corresponding to the central i = 1 . 2 , . . . , m/2. 1, = lm+l-,, frequency of the matching band), with the section impedances (3) varying monotonously from step to step, have been most It can be easily proved that for the structure under considadequately investigated [ 11-[3]. Fig. I(a) shows the structure eration the conditions (3) together with the equations of such a transformer (2, 2 are impedances of the transmission x, = 2, Z2 = x (4) Manuscript received November 5 , 1993; revised February 15, 1996. The authors are with the Central Research Institute of Measuring Equipment, Department of Microwave Structures, 410002 Saratov, Russia. Publisher Item Identifier S 001 8-9480(96)03811-2.

are the necessary and sufficient antimetry conditions. It is generalized in [5] that in order to achieve the global minimum

0018-9480/96$05.00 0 1996 IEEE

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

194

7

I I

z,=z

Z&

z

2

I

cnq---

#I---

---

e, e, (e) Fig. 1. Structures of stepped transformers: (a) Transformer of Class I. (b) Transformer of Class 11. (c) Miniature transformer with the sections of equal lengths. (d) Generalized structure of miniature transformer. (e) New structure of miniature transformer.

of the goal function Ir(0,A)I in the synthesis of any types of with the data obtained in [8] by use of the empirical formula the stepped transformers of cascaded TEM-transmission line sections it is necessary to fulfil the antimetry condition. A0 1 = -arcctg (R The authors [SI have also evaluated the optimum of the 2n solutions obtained. Such an evaluation has been complicated by the fact that the goal function G(A) = ir(Q, A)I is The requirement to fulfil the antimetry condition [SI has multidimensional. With allowance for the device antimetry it is been also confirmed by the results obtained in [6],where the possible to reduce twice the number of independent variables optimum parameters for the stepped transformers of Class I1 of the function G(A), and the vector A dimensions become with maximally-flat (Butterworth’s) characteristics are given. equal to m / 2 . For the simplest case of m = 2 the function Investigation of the m-section transformer of Class I1 ( m = G(A) turns one-dimensional (1-D) (the vector A has only one 2, 4 , . . .) has shown that its amplitude frequency characteristic component AI), which allows us to analyze it both numerically is analogous to that of the m/2-section transformer of Class I and graphically. The numerical analysis of the function G(A1) designed for matching the transmission lines with the same shows that it is multiextremal within the interval (0, 1). It was z and 2 and the same mismatching tolerance /rimax. Yet not found expedient to consider the function G(A1) for A1 > 1 the matching band of Class I1 transformer is only 10-15% since the longitudal dimensions of the transformer turn to be narrower than that of the corresponding Class I transformer, too large under such Al values. Fig. 2 depicts the function and its total length is shorter by a factor of 1.5-2. Besides, G(A1) for Z / z = 2 , x = &/Q1 = rm1.5. As can be seen Class I1 transformers are characterized by a simpler production from the plot, the function G(Al) has two local minimums. technology, due to the fact that only two dimensions of the The first one is located in the interval 0 < AI < 0.2 and is transmission line cross section corresponding to the impedglobal; A1 for such a case is equal to 0.0771, which coincides ances z and 2 are to be realized along their lengths. In case

+ +

mr

195

MESCHANOV et ul.: STEPPED TRANSFORMERS ON TEM-TRANSMISSION LINES

sections ( m is even) of different lengths li and impedances = 1, 2, . . . , m ) . It has been stated in [7] that only structures, for which the relations (7) are true, will have optimum Chebyshev characteristics

z, (i

I\I 0

I

I

I

I

I

0.2

0.4

0.6

0.8

A1

Fig. 2. The goal function of the transformer of Class I1 ( m = 2).

of Class I transformer the number of such dimensions is equal to n 2. Another advantage offered by Class I1 transformers is that it is easy to take into account the effect of electrical nonuniformities arising at the planes of junction of the transmission line sections having different lengths I 1 , 1 2 , . . . ,1, and impedances z and 2 [Fig. l(b)], on the amplitude frequency characteristics. In our case there is only one type of such nonuniformities being caused by jumpwise variation of the geometric dimensions of the cross section of the transmission lines having the impedance z at the regions of their junction with the transmission line sections having the impedances 2. In Class I transformers the number of such nonuniformities 1 [Fig. l(a)], and they are caused by the is equal to TZ jumpwise variation of the geometrical dimensions of the cross sections of the transmission line sections having the impedances 2 1 , z 2 ,. . . , z, at the regions of their junction with the neighbor sections of the transmission line. The transformer length can be further reduced by using the structure described in [9] representing the cascade connection of m transmission line sections ( m is even) of the same length 1 (1 < X,/4), the impedances of which satisfy the following inequalities:

+

+

This transformer structure is shown in Fig. l(c). This transformer is an antimetric device, the antimetry condition being z,z,+1-,

= zZ.

i = 1 , 2 , . . . ,m / 2 .

(6)

The use of this structure allows one to reduce the transformer length by a factor of 2-4 as compared to the analogous device of Class I. The substantial drawback of such a miniature transformer is the necessity to realize a high impedance ratio R, = zmax/zminreaching in a number of cases the values of 30-50. Search of possible ways of eliminating this drawback has led the authors [7] to the generalized structure shown in Fig. l(d). Such a transformer is a cascade of m transmission line

It is easy to prove that the relations (7) are the antimetry conditions for the given structure. This confirms the conclusion made in [5] that the antimetry condition is necessary for achieving the global minimum of the goal function G ( A )= B E ~ ~ 0 2II'(0,A)l 1 during the synthesis of the stepped transformers of all types. The comparison of (7) with the antimetry conditions for the stepped transformers of Class I and IT, as well as with (6) shows that the equations (7) are the generalized antimetry conditions for the stepped transformers of all structure types. Transformer section impedances satisfy the inequality ~,-1

>~

~

> ~. . . -> 231 > Z, > Z,-Z

> . . . > zz (8)

i.e., the impedances of both sections of even and odd numbers decrease in the direction from the transmission line with a higher impedance 2 to z impedance line, the impedance of any section of an odd number being always larger than that of any section of an even number. The following regularity is observed here: the lengths of odd number sections increase in the direction from the transmission line of a smaller impedance z under matching to 2 impedance line, and the lengths of even number sections decrease in the same direction.

111. NEW STRUCTURE FOR MINIATURESTEPPEDTRANSFORMER Another possibility to reduce the stepped transformer length is the use of the stepped structure [see Fig. l(e)], which differs from the one investigated in [ 5 ] ,[6] by its section impedances satisfying the conditions

where z1z2 = xZ,zm< z,xm-l > Z. To solve the optimization problem in form (2), the algorithm based on the linearization method offered by Pshenichni [ lo] has been used. The description of this method as applied to the problem of stepped transformer synthesis is given in the Appendix. In solving the synthesis problem the antimetry property of the given device has been used, which is to be found in case both (9) and (10) conditions are fulfilled

That made it possible to reduce twice the dimensions of the varied parameter vector A . In the general case of msection transformer only one of the impedances ( 2 1 or 2 2 ) and only m/2 - 1 section lengths (e.g., L1, L2, . . . , L,l2-1) should be varied while solving the synthesis problem. The

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

196

z,,n

z2,n

0.1171

52.38

23.86

0.3300

0.0841

72.91

17.14

0.2500

0.0553

114.55

10.91

0.1670

0.0412

155.67

8.03

0.1250

L2,3

0.15

L

0.10

0.05

R

0

TABLE I1 OPTIMUMPARAMETERS FOR THE SIX-SECTION TRANSFORMERS

3%

25

41

0,rad

Fig. 3. Amplitude frequency characteristic of the stepped transformers. Curve I-transformer of Class I; curve 2-new structure of miniature transformer

TABLE I11 OFTIMUM PARAMETERS FOR THE FOUR-SECTION TRANSFORMERS OF EQUAL-LENGTH SECTIONS [9]

lL1,2,3,L1

%(l)

=wzl(2),

Lm/2 =

L/2

'2'*

1 0.065

0.0833

42.38

19.80

63.13

29.49

0.071

0.0625

55.75

13.58

92.03

22.42

0.074

0.0418

82.58

8.45

148.01

12.14

0.1670

1 i.40

0.1250

other impedance value as well as length Lm12 may be defined then from the relations

109.64

0.0313

0.076

I

m12-1

1

'i

I

I

6.17 I

202.48 I

I

L2

-

TABLE IV OFTIMUMPARAMETERS FOR THE FOUR-SECTION

2=1

where L is normalized summary transformer length, which was taken fixed when solving the synthesis problem. Thus the vector of varied parameters is A = ( 2 1 , L I ,L 2 , . . . , Lm/2-1). The results for Chebyshev approximation (2) for the foursection transformer ( m = 4), designed for matching the transmission lines with the impedances z = 25 R and 2 = 50 R in the frequency range of one octave, are given in Table I. Table I1 shows the optimum parameters for the sixsection transformer designed for matching the same lines in the frequency range of one and a half octave. The investigation of the properties of the transformer based on the proposed steplped structure has shown that only the transformers for which the lengths of odd number sections increase in the direction from the transmission line with smaller impedance z to the line with Z impedance, and lengths of sections of even numbers decrease in the same direction, have the optimum characteristics (Tables I, 11). In view of the results given in the present paper, as well as in [5]-[7], we may state that the former regularity of section length variation is inherent in all the stepped transformers designed on the cascade connection of uniform transmission lines of different lengths irrespective of the type of their impedances (whether they assume two alternating values, or are subjected to some other law of variation). Fig. 3 shows the amplitude frequency characteristics for two-section transformer of Class I (curve 1) as compared to those for the corresponding four-section miniature transformer (curve 2 ) both designed for matching the transmission lines with impedances relating as 1 : 2 in the frequency range of one octave; the length of the miniature transformer being 0.125X0, and that of the corresponding Class I transformer

TRANSFORMERS ON GENERALIZED STRUCTURE [7] I

I

I

I

I

I

I

68.73

I

I

0.068

0.0625

0.0833

51.75

18.20

24.15

0.2916

0.070

0.0525

0.0725

62.00

15.28 81.85 20.15

0.2500

0.075

0.0320

0.0510 103.00

10.13 123.39

12.14 0.1660

0.076

0.0205

0.0420 152.90

7.78 160.95

8.18 0.1250

being 0.5X0. For comparison, Tables I11 and IV show the optimum parameters of the known four-section miniature transformers (computed on the data given in [7], [9]), designed for matching the same transmission lines ( R = Z/z= 2) as the proposed transformer (Table I) in the frequency range of one octave. IV. CONCLUSION The investigation of the performances for the stepped transformers of various structures on cascade connection of uniform TEM-transmission line sections has advanced lately. The general patterns of the length and impedance distribution both for the transformers of Class I and Class I1 and those of the generalized structure have been defined. The investigation of the monotonic Class I structures began approximately 40 years ago with the basic works by Cohn, Collin, and Riblet [ 11-[3]. Nonmonotonic structures of Class I1 and the generalized structures have been researched later, primarily by the Russian scientists, with the application of the numerical optimization methods [4]-[7].

MESCHANOV ef al.: STEPPED TRANSFORMERS ON TEM-TRANSMISSION LINES

The use of the new structure ( m = 4) proposed in the present paper allows one to achieve the device summary length equal to 0.125X, (see Table I), which is one-fourth of the length of Class I two-section analog. The same transformer specifications can be provided by the application of the well-known miniature structure [9], offering the same length reduction, but for this it is necessary to realize the following impedances: z1 = 109.64 0 , z z = 6.17 0 , z 3 = 202.49 0, z4 = 11.40 R; Le., the maximum impedance ratio R,,, = z ~ / z z= 32.8. In the structure under consideration the maximum impedance ratio R,,, = 19.39. As compared to the available analogs [ 7 ] , [9], the proposed transformer is the most prospective in terms of productional simplicity, since only four cross-section dimensions corresponding to impedances z , z l , z z , 2 are to be realized along the m-section transformer length. The number of such dimensions in transformers [7], [SI is equal to m 2.

797

problem specificity, i.e., of the fact that the initial problem is that of the continuous minimax. When forming the subscript plurality J 6 ( A k ) , we will take into account only the points, at which the local maximums for a given step are achieved. Then p 1 points will get into the plurality J * ( A k ) at the most, in view of the fact that the number of the local maximums in the nondegeneracy problem (Al) exceeds that of the variables by unity. Problem (A2) is a problem of convex programming. Since its direct solution is difficult, we will transform it into a dual problem. It may be shown that the dual problem is the problem of square programming and has the form

+

II

112

+

under constraints

APPENDIX The linearization method offered by Pshenychni [lo], which is efficient as applied to the problems of discrete minimax, will be used here to solve the problem (2), similar to that described in [5]. First, we will transform (2) into the discrete problem by introducing the set of N >> p points over the interval [ B l , 021 min max F i ( A ) A

l

5

10.26

9 9.375 6 8&

t

9.37

2

u4

f

9.365

' '

9.36

10.18

.........

m710

(double disks) "

'

I

3

"

2

I

"

'

L "

" '

"

'

5

4

10.1 4

6

0

1 o4

1

3

measured

0.08

(a)

1 2

measured

o

0.04 0.06 gap spacing (mm)

0.02

gap spacing (mm) (a)

'OS

0

4

5

6

gap spacing (mm) (b) Fig. 5. (a) Experimental results show a typical resonant frequency tuning of two identical TM modes. The theoretical and experimental results have a good agreement when the gap spacing has a large value. (b) Measured Q-factors of the two TM modes.

C. Coupling Rules and Mechanism

The WG modes are not simply orthogonal eigen-frequency modes. They are hybrid modes and the interaction between the coupled WG modes in the double disk resonators has some interesting characteristics. The coupling between coupled WG modes with the same radial number n = 1 in double disk resonators were experimentally investigated. For isolated double disk resonators, the experiments showed that there are some selection rules for the inductively coupled modes. The inductive coupling between two WG modes occurs only when the differences of the azimuthal and axial mode numbers satisfy the following relations,

a p =O,&l a m = 0, f l where A p and Am are axial and azimuthal number differences of the two WG modes respectively. In isolated double disk resonators, the inductive coupling occurs only between quasiTM modes or between quasi-TE modes. There is no inductive coupling between TM and TE modes. As shown in Fig. 6, the tuning curves with an inductive coupling show a mode transition when the gap spacing increases, i.e., mode TM712

- - - -measured

- - e - - measured

io4 0

0.02

0.04

0.06

0.08

gap spacing (mm) (b) Fig. 6. (a) A typical inductive coupling between two different TM modes (b) Measured 9-factors of the two TM modes.

changes into TM613 and mode TIL1613 becomes TM712 after the interaction. When the two modes are strongly coupled, they become hybrids of each other. The reactive coupling of two modes can occur in a wide frequency range which is larger than their bandwidth. A typical resistive coupling is illustrated in Fig. 7. The coupling only influences the Q-factor of the high Q TM710 mode. This resistive coupling can also occur between TM modes in isolated double disk resonators. Usually resistive coupling occurs inside the bandwidth of the modes. The resonant frequencies of the coupled modes remain unaffected. The above coupling rules are only true for isolated double disk resonators. When the double disk resonators are shielded inside metal cavities, the coupling is modified and become stronger in comparison with the isolated case. In this case, inductive coupling can occur between more sets of modes. For example, the interaction between the TM710 and TE511 in a cavity appears to be a inductive coupling and their coupling coefficient K M 0.0007. However, the modes having close mode numbers, especially the mode numbers which satisfy the relations (9) and (lo), show stronger couplings than the others. The observed phenomenon of the interaction suggests that the couplings between the coupled WG modes are mainly at-

PENG: STUDY OF WHISPERING GALLERY MODES IN DOUBLE DISK SAPPHIRE RESONATORS

9.2

1

9.1 5 n

3

9.1

5. 51

9.05

E

9

LL

L

/

’

853

ACKNOWLEDGMENT The author would like to thank D. Newnnan and A. Gorham for their assistance on machining sapphire disks as well as D. G. Blair, T. Mann, and I. S. Heng for their assistance in preparing this manuscript.

I

TM7,0

8.95

=51,

L 8.gt ’ 0

’

’

’

’

0.05

’

’

’

’

0.1

’

‘

‘

’

0.1 5

gap spacing (mm) (a)

10’

0

0.05 0.1 gap spacing (mm)

0.1 5

(b) Fig. 7. The interaction between two resistively coupled modes. (a) Measured frequencies. (b) Measured &-factors of the two TM modes.

tributed to their electromagnetic field states. This is because the modes with close mode numbers have similar electromagnetic field states that provide large field overlap (or interaction cross section) and thus strong couplings. (This phenomenon may be considered to be analogous to the state transition of atoms in which the transition occurs only between closer states.)

IV. CONCLUSION A mode matching method for determining the resonant frequencies of the WG modes in double disk resonators has been shown to give frequencies accurate to a few tenths of a percent for the fundamental WG modes. The theory gives improved accuracy ( < l % )for calculation of the TM modes with high axial mode number ( p 2 2) in single resonators. Experiments and theory have both shown that the TM modes with a zero or even axial mode number in the double disk resonators have larger frequency tuning than other modes . The interaction between various WG modes in the double disk r e s o n a t o r has been extensively investigated. The coupling of various WG modes is mainly determined by their electromagnetic field states. The closer the mode states, the stronger the coupling of the coupled WG modes. Mode transition may occur between the strongly coupled WG modes.

REFERENCES A. J. Giles, S. K. Jones, D. G. Blair, and M. J. Buckingham, “A high stability microwave oscillator based on a sapphire loaded superconducting cavity,” in Proc. IEEE 43rd Ann. Symp. Freq. Contr., 1989, p. 89. G . J. Dick and J. Saunders, “Measurement and analysis of a microwave oscillator stabilized by a sapphire dielectric ring resonator for ultra-low noise,” IEEE Trans. Ultrason., Ferroelect. Freq. Cont., vol. 37, no. 5 , p. 339, 1990. M. E. Tohar, E. N. Ivanov, R. A. Woode, and J’. H. Searls, “Low noise microwave oscillators based on high-Q temperature stabilized sapphire resonators,” Proc. IEEE Int. Freq. Contr. Symp., 1994. D. G. Blair, E. N. Ivanov, and H. Peng, “Sapphire dielectric resonator transducers,” J. Phys. D: Appl. Phys. vol. 25, p. 1110, 1992. H. Peng, D. G. Blair, and E. N. Ivanov, “An ultra high sensitivity sapphire transducer for vibration measurements,” J. Phys. D: Appl. Phys., vol. 27, pp. 1150-1155, 1994. A. Karp, H. J. Shaw, and D. K. Winslow, “Circuit properties of microwave dielectric resonators,” IEEE Trans. Microwave Theory Tech., vol. MT’I-16, pp. 818-828, 1968. S. Fiedziuszko and A. Jelenski, “Double dielectric resonator,” IEEE Trans. Microwave Theory Tech., vol. MT’I-19, pp. 779-880, 1971. K. A. Zaki and C. Chen, “Coupling of nonaxially symmetric hybrid modes in dielectric resoantors,” IEEE Trans. Microwave Theory Tech., V O ~MTT-35, . p. 1136, 1987. S. W. Chen and K. Zaki, “A novel coupling method for dual-mode dielectric resonators and waveguide filteis,” -l’EEE Trans. Microwave Theory Tech., vol. 38, p. 1885, 1990. M. E. Tobar, “Effects of spurious modes in resonant cavities,” J. Phys. D: Appl Phys., vol. 26, p.- 2022, 1993. G. J. Dick, D. G. Santiago, and R. T. Wang, “Temperature compensated sapphire resonator for ultra-stable oscillator capability at temperature above 77K,” Proc. IEEE International Freq. Contr. Symp., p. 421, 1994. Y. Garault and P. Guillon, “Higher accuracy for the resonance frequencies of dielectric resonators,” Electron. Lett., vol. 12, no. 18, p. 475, 1976. D. Kajfez and P. Guillon, Eds., Dielectric Resonators. Norwood, MA: Artech House, 1986, p. 65. M. E. Tobar and A. G. Mann, “Resonant frequencies of higher order modes in cylindrical anisotropic dielectric resonators,” IEEE. Trans. Microwave Theory Tech., vol. 39, no. 12, p. 2071, 1991. H. Peng and D. G. Blair, “High tuning coefficient whispering gallery modes in a sapphire dielectric resonator transducer,” Proc. IEEE Int. Freq. Contr. Symp., Boston, p. 328, 1994. R. Shelby and I. Fontanella, “The low temperature electrical properties of some anisotropic crystals,” J. Phys. Chem. Solids, vol. 41, pp. 69-74, 1980.

Hong Peng was born in Kunming, China, in 1963. He received the B.S. and M.S. degrees, both in physics, in 1984 and 1987, respectively, from Zhongshan University, Guangzhou, China. He received the Ph.D. degree in physics from the University of Western Australia in 1996. His doctoral dissertation involved the development of ultra-high sensitive sapphire resonator transducers. Since December 1995, he has been with Massachusetts Institute of Technology as a Postdoctoral Associate. Currently his interests are in the detection of dark matter axion, low noise amplifiers, and cryogenic resonators.

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

854

Modeling of General Constitutive Relationships in SCN TLM Leonard0 R. A. X. de Menezes, Member, IEEE, and Wolfgang J. R. Hoefer, Fellow, IEEE

Abstract- The modeling of general constitutive relationships in SCN (symmetrical condensed node) TLM is presented. The technique consists of decoupling the impulse scattering at the nodes from equations describing the medium by using equivalent node sources with state-variable formulation of the constitutive relationships. The procedure requires few modifications of TLM. Numerical examples are presented.

V,'

4

I. INTRODUCTION

I

N THE BASIC TLM formulation, dielectric permittivity and magnetic permeability are modeled by open and short circuited stubs which are connected to the nodes; the characteristic admittance of the stubs is a function of the constitutive parameters [I]. This method is very robust and most appropriate when and are constants. However when the material constitutive parameters are frequency dispersive and nonlinear, the representation by hard-wired stubs becomes computationaly difficult and uneconomical because a modification in time of the stub admittances leads to a modification of the impulse scattering matrix of the node. A better way is to decouple the impulse scattering at the nodes from the equations describing the behavior of the medium by representing the latter by a differential equation or equivalent lumped element network and connecting it to each node by a transmission line of infinitesimal length and characteristic admittance equal to the driving point admittance of the node [2]-[4]. The procedure presented in this paper is an efficient TLM representation of arbitrary constitutive relationships. Using decoupled scattering matrices with equivalent node sources (Thevenin and Norton equivalents of the node), the constitutive equations are expressed in a generic formulation allowing the inclusion of arbitrary medium behavior in TLM models. The resulting equations are solved with the state-variable approach. 11. THEORY This section is divided into four parts. The first describes

the modeling of arbitrary dielectric and magnetic materials in two-dimensional (2-D) shunt and series TLM networks. In the second part, the method is extended to three-dimensional (3-D) SCN TLM. In the third part, the networks describing Manuscript received June 10, 1995; revised February 15, 1996. This work was supported in part by the Brazilian Govemment agency Conselho Nacional de Pesquisa (CNPq). The authors are with NSERCMPR Teltech Research Chair in RF Engineering, Department of Electrical and Computer Engineering, University of Victoria, Victoria B.C., Canada. Publisher ltem Identifier S 0018-9480(96)03804-5.

Circuit

Fig. 1. Two-dimensional TLM shunt node NI

Fig. 2. Two-dimensional TLM shunt node N2.

the medium behavior are formulated in terms of equivalent node sources, and in the fourth part the resulting differential equations are expressed and solved using the state-space approach. A. Two-Dimensional TLM Formulation

General isotropic materials can be modeled in 2-D-TLM by reactively loading each node of the network. As mentioned above, the reactive load can be either modeled by a reactive stub (Fig. 1) or by a reactive lumped element network (Fig. 2). The latter formulation will now be given for general nondispersive linear media, for the sake of simplicity. Since the characteristics of shunt and series-connected TLM networks are related to each other by duality, we will only

0018-9480/96$05.00 0 1996 IEEE

DE MENEZES AND HOEFER: MODELING OF GENERAL CONSTITUTIVE RELATIONSHIPS IN SCN TLM

derive the formulation for the shunt case and simply state the analogous results for the series case. In the shunt connected TLM cell the reactive network across the node models the polarization of the medium in the presence of an electric field, and the current flowing through it is the polarization current. The equations relating this current to the electric field describe the dielectric response of the medium. Maxwell's second curl equation in two dimensions ( d / d g = 0)

855

it with the scattering mechanism in the TLM network, the total voltage and current across the reactance must be expressed in terms of incident and reflected voltage impulses. To this end, the shunt reactance is connected to the node via a transmission line of infinitesimal length and characteristic admittance Y,. In order to avoid multiple reflections on this line, Y, is matched to the driving point admittance of the node, i.e. Y, = 4Y0, where YOis characteristic admittance of the link lines. The voltage and current across the shunt reactance can thus be related to the voltage impulses incident and reflected at the node on this transmission line as follows:

is modeled by the shunt type 2-D-TLM network as

+ wy:

v, = v:,

The scattering matrix of the node is

with the equivalences at

=-Hx

EO

=2co

2,

E,

-H,

v,

2co

G

E,

+ c,

where COis the capacitance per unit length of the link lines. For a square cell with Ax = Az = Althe total displacement current is thus dv av ad = i d 0 f i, = 2cOa12 f c,ai-! (3)

at

and the reflected voltage v,T, is calculated as

at

where Zd is the total displacement current, ido is the displacement current in vacuo and i, is the polarization current in the medium. The equivalence between

v:, =

i(v;

+ w; + w; + v l ) .

used in (8) is calculated Therefore, the incident voltage using (9) with (7) and (6). In the series node case, the decoupled scattering matrix is 'u1

3 1

€0

-1 0

I

= 2Co c,= 2cox = 2Co(t, - 1)

+ i o

0

obtaining for the polarization current in (3) i, = 2CO(€,- 1 ) AdV lL.

(4)

at Furthermore

Co =

gm=; Yo

=

yoat

~

a1

1 1 - 1 - 1 1

;' ; ;]

3

[;$=:[:l yields with

0

0

0 .0 0 1 0 0 0 -1.-

l o

The current and voltage over the node a.re defined as

vmag(t)= v & ( t ) (5)

where LO and YOare the inductance per unit length and the characteristic admittance of the link lines, respectively, At is the timestep, and c the speed of light (on the link lines). Hence

where p , is the normalized polarization P y EE P y / E O .

The shunt reactance connected to the node is usually described by a current-voltage relationship. In order to connect

(9)

i(t) =

.k(t)

+ &(t) -

vk(t)

420

(1 1)

and the magnetization vector is

where (in the linear nondispersive case) m ( t ) = ( p T l)i(t) is the normalized magnetization vector, wrnag(t) is the magnetization voltage and V; ( t ) v: , ( t ) are the reflected and incident voltages at the input port of the network modeling the magnetic behavior of the medium.

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

856

.;

[ 6 ] ,with the same numbering scheme used in [7]

B. Obtention of the SCN Formulation The SCN formulation is obtained by combining 2-D shunt and series nodes [6]. These nodes are decoupled from the medium using the procedure described in Section 11-A. The resulting medium constitutive equations are

%

YOAl ipolz

= (2T)

A1

ZoAl

umagx =

( 2 7 )

ZoAl

dt

=( 2 T )

vmagy

my(C

260

(a+YOAl

=Pz(t)g = Yr(u:z

-

= Yr(&

-

vk,)

= yr(u:z

-

= ( 2 7 )

dm,

dt

Mr(t)-2

uz = u:, i

i,(t) = (wi(t) - u k ( t ) - $(t)

,i z )

-

-

The conditions (15) to (16) can be expressed in the system

,

v, =

- vi@)) - &,(t))/4

-

+ ukz

4,)vy = viy + 4,

.az,

+ w$(t)

+ v i ( t ) w&(t)) - wLy(t))/4 i z ( t )= ( ~ i ( t )~ i ( t+)vi,(t) - ~ ; l ( t ) )~ & , ( t ) ) / 4 . (16)

i y ( t )= (w;(t)

(13)

iy iz)

Py = p z (wz,,uy 7 w.2, i z , 2 , ipolr

+ W ; ( t ) + 'Ui(t)+ &(t))/4 f v;,(t) u y ( t ) = ( W i ( t ) + wj(t) + + 'U;l(t))/4 + vby(t) v*(t)= (vk(t) + ? J ( t+) w $ ( t ) + v4,(t))/4 + &(t)

w,(t) = ( u f ( t )

Vg(t)

A1

mz(t) =

p , = p , (vz, wy 1 u z , i, ipoly

with Y,/Yo = Zr/Z0= 4 and

dm,

dt

2

ZoAl

dPd

A1

ipolz

= ( 2 T )

= My(+

Vmagz

P Z ( t )

+ +

(15)

A1

A1 P,(t) = P z ( t ) =

dm,

dt

A1 2

YOAl d P y

ipolz

w; = wZ(t) + z,iz(t)- v$(t) u; = w,(t) - Z r i y ( t )- w; w; = vz(t) Zriy(t) - v; u; = u r ( t ) - Z&(t) - uto wy0 = W Z ( t ) Z r i y ( t )- 7J; 21; = u y ( t ) + Z,i,(t) w; = w z ( t ) - Z r i z ( t )- wg ?JYl = u y ( t ) - Z,i,(t) - w; w; = wy(t) - Z&(t) - w; u8' = wy(t) + Z r i z ( t )- w; wT2 = wz(t) + Z,iz(t) - w: -

mz ( t )= M z ( t )-

Pz ( t ) = p, ( t ) ipoly

= u, ( t )- z,i,(t)

[UT]

=

+

[S][Vi] i[w$] [vi]= [T][vi] (17)

where, as shown in (18) and (19) at the bottom of the page, and

+ 4,

and the subscript t denotes the transpose vector. The reflected voltages w; and w: are calculated as follows:

+ v; + w; + w&) &, = (vi + w; - w8) + w; + + &) ?JLy= (vi w; + wi,) u,Ty = = ;(w; + v; + w 7 2 + vio) wkr = (wf - w; - w;, + w;,,.

we', =

1

- ?Jk

1

7J;

the reflected voltages on branches 1-12 are calculated using

-0 1 1 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 1 0 0 0 [SI= 1 4 0 0 0 - 1 0 0 1 - 2 1 - 2 0 0 0 - 1 0 0 - 1 0 - 2 1 _ - 2 1 - 1 0 - 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 [TI = 1 4 0 0 0 1 - 1 0 - 1 0 0 0 -1 0 - 1 0 0

0 0 0 1 0 1 2 1 0 1 0 0 0 0 1 0 1 0

-

?J; -

(21)

0 0 1 0 - 1 - 2 0 0 - 2 - 1 0 1 0 1 0 0 - 2 - 1 0 - 1 - 2 0 0 1 0 1 - 2 - 1 0 1 0 0 0 1 0 - 1 - 2 0 0 1 0 1 0 1 0 0 0 1 0 0 0 1 0 - 1 0 0 0 1 0 1 - 2 1 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 1 0 0 0 1 0 0 1 0 1 0 0 1 0 1 0 0 1 0 0 1 - 1 0 0 0 0 0 0 1 - 1 0 0 0 0 0 0 1-10

1 0

DE MENEZES AND HOEFER: MODELING OF GENERAL CONSTITUTIVE RELATI(3NSHIPS IN SCN TLM

857

2Al d d dt

is,@)=22,~k,(t) = 8vk,(t) = - -m,(t)

+ 4i,(t). (25)

i

The sources are calculated at each timastep by

+ + W i ( t ) + v";(t)) + + + vi&)) isr(t) = 2(v5(t)+ vk(t) + v i ( t )+ & ( t ) )

i s z ( t )= 2 ( v l ( t ) vk(t)

i s y ( t )= 2 ( 4 ( t ) v%(t) v i ( t ) (a)

(b)

us,@) = 8(vS(t) - v$(t) + v;(t) - vi@))

Fig. 3. Norton equivalent source of the node (dielectric behavior).

Zisy(t)

The link between the incident and reflected vectors given in (20) is obtained by (13) and (14).

C. The Equivalent Node Sources The disadvantage of the substitution of (14) into (13) and subsequent discretization of the resulting expression is the loss of flexibility of the constitutive relationships modeling, since in this approach the final discretized expression will be different according to the medium equations. Therefore, it is important to deduce a robust solution procedure that enables (13) to be used independently of the constitutive equation. This is done considering the Norton and Thevenin circuit equivalent of the node as sources connected to the networks. The voltage and current on the input port of each network can be expressed as 1) Shunt network:

ipo,(t) = 2Yrv,.(t) - Y r v ( t ) v i ( t ) = v ( t )- v I ( t ) . (22) 2) Series network: = 2z,w;(t)-Z,i(t)

v i ( t ) = -(4i(t)-vK(t)). (23) These equations can be rearranged and substituted into (13) resulting in ?Ima,&)

A1 d -p(t) + Yrv(t) c dt A1 d us@)= 2ZrvK(t) = 2 0 2 -m(t) Zri(t). (24) c dt The representation of (24) in the equivalent circuit form for the polarization is shown in Fig. 3. Applying the same procedure to all directions, the equivalent current sources are obtained from the polarization current i p o l and the total voltage across the node v ( t ) ,while the equivalent voltage sources are deduced from the magnetization voltage wmag and total current over the node i ( t ) i s ( t )= 2YrvL(t)= Y02-

+

2Al d i s x ( t )=2Y,~,',(t)= 8v,',(t) = - - p x ( t )

+ 4v,(t)

isy(t)

+ 4vy(t)

d dt 2Al d =2Y,v;,(t) = 8 ~ i , ( t )= - - ~ , ( t ) d dt

= 8(v;(t) - v i ( t )

+ v i ( t )- vio(t)) +

vsZ(t) = 8 ( v i ( t )- v i ( t ) - vj2(t) vil@)).

(26)

The solution of (26), (25), and (15) will result in the reflected voltages in lines 1-12. The propagation between nodes is not affected. The incident voltages (20) can be calculated using the results from (25) substituted into (14). This procedure is applicable to all kinds of constitutive relationships. The adaptation to usual TLM programs, [7], is done by setting Y, = Y, = Y, = Z, = 2, = Z, = 4, obtaining the reflected voltages for the stubs, using (25) and (14) to obtain the incident voltages from the stubs, and then calculating the reflected voltages in branches 1-12 with the scattering matrix

(18). D. The State-Variable Approach The use of equivalent node sources allows the solution of the network using nodal or Tableau analysis, [9]. Therefore, a SPICE circuit simulator could be used to solve the network at all nodes at each timestep. The problem of this approach is the need to formulate an equivalent circuit of the medium. Although this may be an easy task for most linear dispersive isotropic materials, that is certainly not the case for an arbitrary constitutive relationship. However, the state-variable formulation of the constitutive relationship equations avoids this problem. In this approach, the use of equivalent circuits to model the constitutive relationships is not discarded but it is not necessary, since the state-equations can be obtained directly from nonlinear differential equations. The state-variable technique is easily incorporated into the TLM simulator without loss of generality. A general procedure for linear differential equations shown in [5] is outlined for a medium described by a fourth-order linear differential equation a-d4f +b- d3f + e -d2f + d - df -+ef = g

dt4

dt3

resulting in

where z l ( t ) = f ( t ) .

dt2

dt

(27)

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

858

The resulting equations can be solved either analytically or numerically. The analytical solution of (28) would result in a discretized convolution procedure, restricting the use of the formulation to linear materials, since the convolution procedure requires the linearity of the system. If the numerical approach is used, several discretization schemes may be chosen. However, there are two major schemes that are very attractive, for reasons of stability and efficiency: 1) Backward Euler scheme [lo] 2) Approximate Trapezoidal scheme (or first order Pad6 approximant) The first scheme introduces losses in the final result, but there is no frequency shift, and it is very simple to implement. The second scheme is far more precise than the former, conserves energy and uses the same kind of discretization as used in TLM, Therefore, for linear isotropic nondispersive materials, the results given by this approach and usual TLM are virtually indistinguishable. Both schemes transform the continuous state equations d [74t)l dt [z(t)I = [AIMt)l + [BI [Y(t)l = [CI[4t)l + [DlMt)l

The state-variable form is

(34) The final equations are

where [U] is the 3 x 3 unitary matrix, [w(t)] and [i(t)] are the vectors containing the voltages and currents in z , y and z ; [wi(t)] and [ w L ( t ) ] are the reflected voltages given by (20) and

B. Linear isotropic dispersive medium

-

(29)

into the discretized form

In the case of a dispersive dielectric medium modeled by a first-order Debye approximation, the frequency domain permittivity function is [ 111

[z(t+ at11 = [PI[z:(t)l + [QI iu(t)I + [RI[ x ( t - At)] [y(t At)] = [C][ x ( t + a t ) ] + [D][u(t+ At)]. (30)

+

The matrices [PI,[Q], and [R]will depend on the discretization scheme 1) Backward Euler scheme:

(37) The relationship between E and D in this dispersive material will be modeled by the RC circuit shown in Fig. 4, with the analogies c 1 = 2At(E,

[PI = ([U] - At[A])-l

[&I

=

(VI

-

R=

WAl)-'At[Bl

[RI = [OI

(31)

[Ql

=([UI

at

+WAI)

WAI)-'T[E a t [R] = ([U] - At[A])-'--[B]. 2

1)

c,= 2&(€,

- E,)

2At(€,-.),t

The state equation describing the circuit will be

where [U] is the n x n identity matrix, and [0] is the null matrix. 2) Approximate trapezoidal scheme:

[PI = (PI- at[Al)-'([Ul

-

70

d = [AI[%] -[vu] dt [4u1

+ [BI[.UeuI + Pl[4,1

=[CI[4

(38)

with

-

(32)

111. FORMULATION OF THE STATE-VARIABLE EQUATIONS This section shows the implementation of the state-variable equations for several kinds of media. The presented results are the continuous state-variable equations.

where u = x , y and z , uu,,,(t) is an auxiliary variable used in the state-equation description of the system, Y,is equal to 4 and viu,viu are the reflected and incident voltages at the A. Linear isotropic nondispersive medium input port of the network in the u direction. Therefore for each In this case electric field component ( x , y and z ) , a system of equations hl d . (39) has to be solved. ~ s ( t =) 2 Z r ~ & ( t = ) Zo(br - 1)2- - ~ ( t ) Z r i ( t ) The dispersion analysis of TLM using state-variable apc dt A1 d proach shows that the timestep should be at least A t < r O / l O O i s ( t )= 2Yrv3t) = YO(€,- l)2-v(t) Yrv(t) ( 3 3 ) c dt to obtain accurate results. If a small frequency shift is allowed

+

+

~

DE MENEZES AND HOEFER: MODELING OF GENERAL CONSTITUTIVE RELATIONSHIPS IN SCN TLM

the limit can be decreased to At < ~ / 2 0 .This restriction applies to backward Euler and approximate trapezoidal discretization. In the case of second order materials, [12], the restriction is the same for the approximate trapezoidal case and it is worse for the backward Euler case. This rule of thumb for the discretization is valid as long as E , , p, < 20.

Considering the anisotropic material with nondiagonal tensor -

--

u

Fig. 4. Equivalent circuit model of a first-order Debye dielectric.

C. Linear anisotropic nondispersive material

P = [FIE = ([E]

859

Comparison of Results - WR28 waveguide

[U])E

= [GIH = ([PI -

(40)

[WH

where Exx

[:I

[U]= 0

Exy

Pxx

Ex,

Pxy

Pxz

:]

1 0

.

frequency (GHz)

The constitutive state-equation is

d 2Atdt

This relationship is obtained directly from (40) and (41) without the need for an equivalent circuit.

D. Nonlinear material Consider the constitutive relationship of a nonlinear medium

+

~ ( t= )( e , - 1 ) ~ ( t ) T ( E ( ~ ) ) ~ (43) This is an example of a self-focusing material, because the effective dielectric constant increases with the amplitude of the wave. An electric field propagating in a waveguide loaded with a dielectric strip with this constitutive relationship tends to be more concentrated in the strip with increasing field values. Above certain power levels one observes the formation of spatial solitons in the guide. The state-equation will be d Yr -v(t) = v(t) dt %[(E, - 1)2Al/c 3 ~ ( v ( t ) ) ~ l

+

Fig. 5. The cutoff frequency of a WR-28 waveguide filled with dielectric. The frequency is obtained after a Fourier transform of the time domain response of an impulsive excitation of a cavity with the dimension of the guide. Results obtained with three different methods are compared. (a) stub-loaded SCN-TLM. (b) State-variable equations (backward Euler discretization). (c) State-variable equations (approximate trapezoidal discretization).

1) Comparison between cutoff frequency results obtained for a dielectric-filled isotropic waveguide using statevariable TLM and stub-loaded TLML The example was a WR-28 waveguide (7.1 12 mm by 3.556 mm, backed by magnetic walls) with regular mesh with a discretization of 24 x 12 x 4 filled with a dielectric with E , of 2.22. The first dominant mode of the cavity has the same frequency as the WR-28 guide. The results are shown in Fig. 5. 2 ) The calculation of the scattering parameters of a parallel plate waveguide with an air/dispersive material junction. The parallel plate waveguide was modeled by a mesh of 200 x 10 x 5 nodes (14.65 x 0.7325 x 0.36625 mm) with the dielectric constant of E, = 1 as shown in Fig. 6. In the first case (Fig. 7), the dispersive dielectric was modeled as a first order Debye medium with parameters seconds. In E, = 1 . 8 , ~ =~10.0 and TO = 9.4 x the second case (Fig. 8), the dielectric was modeled as a second order Debye medium with t, = 2 0 . 0 , ~=~ 60.0, fo = 5 GHz and 6 = 0.3. Both results were calculated using state-space equations discretized using backward Euler scheme. 3) The calculation of the cutoff frequencies of a sapphire filled WR-28 waveguide with the same discretization used in the first example. The permittivity tensor is ( t ucos'

cp

+

tw sin2 c p )

sin 2cp

;(eu

--

(ewcos2 (p

t,)

"I

sin 2cp sin2 'p)

+ E, 0

€2

(45)

IV. NUMERICAL RESULTS This technique was validated by comparing SCN TLM results for several materials:

where the dielectric was sapphire 11.49), [8].

(E,

= 9.34

c,

=

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

860

-

Analytical Cutoff Frequency

Axis Angle

SCN TLM

Error (%)

6.2221 GHz

O0

6.21 GHz

0.19

6.5354 GHz

450

6.57 GHz

0.53

6.9012 GHz

90’

6.90 GHz

0.02

Reflection Coefficient

Electrk Wafl Absorbing Wall Electric Wall

DEreceion of P Fig. 6

~

~

~

~

~

~

i

Parallel plate waveguide half-filled with dispersive dielectnc

#Fig. 8.~ Scattering parameters from an aiddispersive dielectric transition. Solid-line exact result; dashed-line: result calculated with TLM. The RMS error over the whole frequency band is in the range of 3.0% (with a maximum of 6.0% at 36 GHz).

Reflection Coefficient

O

82x exact solutions for the anisotropic case. Good agreement was observed in both cases.

REFERENCES [ l j W. J. R. Hoefer, “The transmission line matrix method-theory

[2] [3]

frequency (GHz) [4j Fig. 7. Scattering parameters from an air/dispersive dielectric transition. Solid-line: exact result; dashed-line: result calculated with TLM. The maximum error over the whole frequency band is in the range of 0.4%.

[5] [6]

The optical axis lies on the xy plane and was rotated by an angle cp with respect to the 1c axis. The problem was calculated for cp of 0”,45”, and 90” using backward Euler discretization. The comparison between the exact and calculated results is shown in Table I.

[7] [8j [9]

V. CONCLUSION The technique presented in this paper can be used for modeling general constitutive relationships and requires few modifications to a TLM program. A general description of the medium relationships was obtained with equivalent node sources and the state-variable approach. The technique was validated by comparison with stub-loaded SCN results and

[lo]

[ 111

[ 121

and applications,” IEEE Trans. Microwave Theory Tech., vol. MTT-33, pp. 882-893, Oct. 1985. P. Russer, P. P. M. So, and W. J. R. Hoefer, “Modeling of nonlinear active regions in TLM,” IEEE Microwave Guided Wave Lett., vol. 1, pp. 10-13, Jan. 1991. L. de Menezes and W. J. R. Hoefer, “Modeling frequency dependent dielectrics in TLM,” in IEEE Antennas Propagat.-S Symp. Dig., June 1994, pp. 1140-1143. ~, “Modeling nonlinear dispersive media in 2-D-TLM,” in 24th European Microwave Conference Proc., Sept. 1994, pp. 1739-1744. R. DeCarlo, Linear Systems-A State-Variable Approach with Numerical Implementation. Englewood Cliffs, NJ: Prentice-Hall, 1989. P. Naylor and R. Ait-Sadi, “Simple method for determining 3-D TLM nodal scattering in nonscalar problems,” Electronic Lett., vol. 28, pp. 2353-2354, Dec. 3, 1992. P. B. Johns, “A Symmetrical Condensed Node for the TLM Method,” IEEE Trans. Microwave Theory Tech., vol. MTT-35, n. 4, pp. 370-377, Apr. 1987. N. G. Alexoupoulos, “Integrated-circuits structures on anisotropic substrates,” IEEE Trans. Microwave Theory Tech., vol. MTT-33, no. 10, pp. 847-881, Oct. 1985. L. Chua, C. Desoer, and E. Kuh, Linear and Nonlinear Circuits. New York: McGraw-Hill, 1987. W. Press, B. Flannery, S. Teukolsky, and W. Vetterling, Numerical Recipes in PASCAL. Cambridge, MA: Cambridge Univ. Press, 1990, ch. 4. R. Luebbers, F. P. Hunsberger, K. Kunz, R. Standler, and M. Schneider, “A frequency-dependent finite-difference time domain formulation for dispersive materials,” IEEE Trans. Eletromagn. Compat., vol. 32, pp. 222-227, Aug. 1990. R. Luebbers and F. Hunsberger, “FDTD for nth-order dispersive media,” IEEE Trans. Antennas Propagat., vol. 40, pp. 1297-1301, Nov. 1992.

DE MENEZES AND HOEFER: MODELING OF GENERAL CONSTITUTIVE RELATIONSHIPS IN SCN TLM

Leonard0 R. A. X. de Menezes (S’89-M’91) was bom in Brasilia, Brazil, in 1966 He received the degree in electrical engineenng in 1990 and the M.Sc. degree in 1993, both from the University of Brasilia (UnB) He is the recipient of a scholarship of the Govemment of Brazil, and is currently working toward the Ph.D. degree in electrical engineenng at the University of Victoria, Victoria B C., Canada His research involves TLM modeling of general media, with emphasis on frequency-dependent and nonlinear dielectric material. His interests are in the areas of numerical modeling, radio frequency components and systems, semiconductor devices, and electromagnetic theory.

86 1

Wolfgang J. R. Hoefer (€7’91)received the Dip1.Ing. degree in electncal engineering from the Technische Hochschule Aachen, Germany, in 1965, and the D Ing. degree from the University of Grenoble, France, in 1968 During the academic year 1968 and 1969 he was a Lecturer at the Institut Universitaire de Technologie de Grenoble and ai Research Fellow at the Institut National Polytechnique de Grenoble, France In 1969, he joined the Department of Electrical Electncal Engineenng, the University of Ottawa, Canada where he was a Professor until March 1992. Since Apnl 1992 he has been the NSERCMPR Tilted Industrial Research Chair in RF Engineering in the Department of Electrical and Computer Enpmeenng, the University of Victoria, Canada Dunng sabbatical leaves he spent six months with the Space Division of AEG-Telefunken in Backnang, Germany (now ATN), and six months with the Electromagnetics Laboratory of the Institut National Polytechnique de Grenoble, France, in 1976 and 1977. In 1984 and 1985 he was a Visiting Scientist at the Space Electronics Directorate of the Communications Research Centre in awa, Canada. He spent a third sabbatical year in 1990 and 1991 as a V ng Professor at the Universities of Rome “Tor Vergata” in Italy, Nice-Sophia Antipolisi in France, and Munich (TUM) in Germany His research interests include numerical techniques for modeling electromagnetic fields and waves, computer aided design of microwave and millimeter wave circuits, microwave measurement techniques, and engineering education. He is the co-founder and Managing Editor of the International Journal of Numerical Modeling. Dr Hoefer is a Fellow of the Advanced Systems Institute of British Columbia

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

862

High-Frequency Reciprocity Based Circuit Model for the Incidence of Electromagnetic Waves on General Circuits in Layered Media Frank Olyslager, Member, IEEE

Abstract- Traditionally a circuit on a high-speed multichip module (MM) or a microwave monolithic integrated circuit (MMIC) is represented in an equivalent circuit by S-parameters for the different components, such as filters or bends, and by transmission lines for the interconnections between the components. Nowadays the S-parameters of the components are easily determined by a numerical electromagnetic analysis. Different components close to each other will interact, often this interaction is unwanted. In the present contribution we develop a circuit model for these interactions without having to perform a global electromagnetic analysis of the interacting components. These interactions are then represented by discrete and distributed sources in the equivalent circuit. Our technique is based on reciprocity and is focused on the surface wave interaction which is often the most important one. Each component is characterized by a surface wave radiation pattern.

I. INTRODUCTION ULTICHIP modules (MCM’s) and microwave monolithic integrated circuits (MMIC’s) typically consist of metalization patterns embedded in a layered structure. Certainly for a MCM, and often also for a MMIC, it is not possible to perform a global electromagnetic simulation that incorporates all the interactions between the different circuits on a MCM or MMIC. Modeling large parts of a circuit at once is very CPU-time consuming and often impractical because every time the design of the circuit is changed the modeling has to be repeated. Due to the layered nature of the substrate the most important interaction between different separated circuits on the same MCM or MMIC is surface wave coupling. Our aim is to characterize each circuit by a surface wave radiation pattern and to represent the interaction of a surface wave with a circuit by a current and voltage source in the circuit model of the circuit. In essence we approach the problem as a transmitter and receiver surface wave antenna problem. A typical circuit on a MCM or MMIC consists of interconnections and what we will call components. These components are everything which deviates from an interconnection structure such as bends, steps in width, filters, active components, lumped elements, via holes, air bridges, etc.. In a high frequency circuit description the interconnections, which act as waveguides, are represented by an equivalent Manuscript received June 12, 1995; revised February 15, 1996. The author is with the Department of Information Technology, University of Ghent, 9000 Ghent, Belgium. Publisher Item Identifier S 0018-9480(96)03803-3.

transmission line model and the components are represented by their S-parameters. In the past much theoretical effort has been spent in constructing equivalent transmission line models for high-frequency interconnections. In [ 11 a rigorous equivalent transmission line model has been derived based on reciprocity considerations. In the same publication the meaning of the impedance level of the transmission line model and at the same time the meaning of S-parameters for connected components has been carefully investigated. In [2], based on the same equivalent transmission line model and on the Lorentz reciprocity relation, a circuit model for the incidence of electromagnetic waves on interconnections has been investigated. The impinging wave is represented as distributed current and voltage sources in the transmission line. Traditionally, components are described by their S-parameters at the ports where the electromagnetic field is assumed to be modal. Nowadays electromagnetic field simulators are available to determine the S-parameters of a very large variety of components. Although we will mainly concentrate on surface wave coupling, it will be shown that other coupling mechanisms such as space wave coupling can be handled with the same techniques. Because of their importance and to make the theory more intelligible we will start with planar perfectly conducting circuits. Later the generalization to general three dimensional, not necessary perfectly conducting, objects embedded in layered media will be discussed. This makes the theory also applicable to dielectric waveguide circuits. 11. SCHEMATIC REPRESENTATION

To focus our attention let us consider the planar metalization structure of Fig. 1. The figure shows the top view of the structure and it is assumed that the metalization is located inside or on top of a stack of layers which can be backed with a ground plane. The structure consists of two separated circuits, indicated with A and B. Each circuit consists of a number (two in the case of the figure) of interconnections (hatched regions on the figure) and a component part. The interconnections, indicated by the subscripts ‘ I ,l’,‘ I ,2 ’ , etc., provide interaction with the external world, Le., with other circuits. In the sequel we will assume two interconnections for each component but the theory is of course valid for any number of interconnections. The component part, indicated by the subscript ‘C’ is an irregular planar metalization pattern,

0018-9480/96$05.00 0 1996 IEEE

OLYSLAGER: HIGH-FREQUENCY RECIPROCITY BASED CIRCUIT MODEL

863

Fig. 1. Interaction between the component parts of two planar metalization structures and the corresponding representation by discrete sources.

such as a filter, a step in width or just a bend. The places where the interconnections are connected to the component part are called the ports of the component. Suppose for a moment that only circuit A is present and that we have modeled this structure with an electromagnetic simulator. From this simulation we know the current densities on the metalization for each excitation of the structure by a fundamental mode on the interconnections AI,^ and AI,^. It is assumed that the fields along the interconnection parts AI,^ and AI,^ are monomodal and that other higher order modes, discrete or continuous, generated in the component part A c have died out and are negligible compared to the fundamental modes at the ports. If this is not the case then the interconnection part must be reduced and the component part must be enlarged. From the electromagnetic simulation we obtain a circuit model for circuit A consisting of two transmission lines, which represent the interconnections, and a scattering parameter matrix S A which describes the component part A c . It has to be remarked that if the interconnection propagates more than one fundamental mode (for example an even and odd mode in two parallel microstrips) or if also higher order modes are important, for example at higher frequencies, the interconnections can be represented by a coupled set of transmission lines. The generalization of the theory to multimode interconnections is presented in Appendix A. A similar electromagnetic simulation for the circuit B

yields the surface current densities on the metalization and the scattering parameter matrix S B . If both circuits are present they will interact. This means that the whole structure starts to act i i S a four-port. One way to characterize this structure is to perform one global electromagnetic simulation for the whole structure and derive the 4 x 4 scattering parameter matrix of the structure. This might still be manageable for the simple structure depicted in Fig. 1 but this technique is certainly to CPU-time consuming for a whole MCM. To overcome this problem we will adopt a different technique to take into account the interaction between circuit A and B. However, we have to rnake the assumption that the effect of the presence of circuit B on circuit A is negligible when circuit B is not excited at its interconnections. Clearly this assumption will fail, e.g., if circuit A and B are parts of the same filter structure which is designed to have strong interaction between its different parts. However, when circuit A and B are truly two different structures then the present method will be able to calculate the electromagnetic interference between both circuits. For planar structures in layered media the most important interaction between such different structures is due to surface waves. These surface waves are eigenmodes of the layered stnrcture in the absence of the metalization patterns and they are excited by the current densities on discontinuities, such as the components A c and Bc. Except for this surface wave interaction, there is also

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

864

e e

BIJ

73

-

-

lB,l(s)ds

'B,{

SB

interaction through space waves, Le., through radiation in the air above the metalization. In this case Ac and Bc act as real antennas. However, it is often the surface wave interaction which is dominant because the surface waves decay more slowly than the space waves as a function of the radial distance. Now we will give a schematic overview and circuit representation of the interaction between circuits A and B. We will separate the interaction between structure A and B into four parts. The first type of interaction together with the circuit model is shown in Fig. 1. This first type of interaction takes into account the incidence of electromagnetic waves, generated by the currents of A c , on Bc and vice versa. So if circuit A is excited by a mode along one of its interconnections, the currents on A c will generate surface waves. These surface waves will induce currents on Bc and excite the fundamental eigenmodes, which propagate away from Bc,in the interconnections B I , and ~ B I , ~In. the circuit model this interaction is represented by the current and voltage , , ~VB,,, , and at the ports of structure B sources I B , ~I B which depend on the currents and voltages i ~ , 1i ,~ , 2V, A , ~ , and VA,2 at the ports of circuit A. Indeed Z A , J , ~ A , ~ , V A , J and v A , ~determine or represent by what amount the eigenmodes in the interconnections of structure A are excited

v,,~

e

BC \B,2

1B,2(s)ds

-

B1,2

e t-*

and I B ) I~ B, ; ~ V ,B , ~and , V B ,represent ~ the excitation of the eigenmodes in the ports of structure B. The sources I B , and ~ VB.~(IB,~ and V B , ~are ) designed such that they generate a wave or mode in the transmission line B1,1(B1,2)which propagates away from the component Bc. Due to linearity I B, , v~B,, 1 , and V B ,are ~ linear combinations the sources I B , ~ of the circuit parameters Z A , J , i ~ , 2V, A , J , and W A , ~The . above reasoning can be repeated when surface waves generated in B c interact with A c . Fig. 2 depicts the second type of interaction. The surface waves generated at A c will not only generate eigenmodes in the interconnections B I , and ~ B I ,through ~ interaction with the component part Bc of circuit B but also directly through in~ B I , themselves. ~ teraction with the interconnections B I , and The interaction of impinging waves on interconnections was studied in detail in [2].The circuit equivalent of this type of interaction is also shown in Fig. 2. The interaction is represented by distributed sources I ~ , l ( sd)s , I B , z ( ~d)s , v B , l ( s ) ds, and VB,~(S) ds in the transmission lines which depend on the currents and voltages ~ A , Jz A, , 2 , V A , J , and w A , 2 at the ports of A c . Here s denotes the distance from the component along the considered transmission line. Of course by a simple integration it is possible to concentrate the distributed sources along a part of a transmission line in a discrete source placed at some

OLYSLAGER: HIGH-FREQUENCY RECIPROCITY BASED CIRCUIT MODEL

iA, 1

e

865

vA,l -1

A, SA 0

0

e--,

U

ds

ds iB,dS)

vB, 1

‘B,2

-1

L

iB,2(S)

BC

SB 0

0 U

ds

e-,

ds

Fig. 3 . Interaction between the fields generated at the interconnection parts with the component part of the other structure and representation by discrete sources.

point of that part. This discrete source will have the same effect outside the considered part of the transmission line as the distributed sources on that part. Normally an eigenmode propagating along an infinite interconnection or waveguide does not interact with surface waves. This is however not true for a finite or semi-infinite interconnection line. This brings us to a third type of interaction which is the reciprocal of the previous type. The interconnections AI,^ and AI,^ will generate surface waves which, through interaction with Bc will generate eigenmodes and B I , ~This . situation is shown in the interconnections B I , ~ in Fig. 3 together with the circuit equivalent. The circuit equivalent consists of discrete sources I B , ~I B , ,V~B,, ~ and , V B , at ~ the ports of Bc which depend on the currents and voltages i A , J ( s ) , ~ A , z ( s WA,J(S), ), and V A , ~ ( S ) of each elementary part ds of the interconnections AI,^ and AI,^. Of course one could again integrate along a part of the interconnections AI,^ or AI,^ to concentrate the effect of this part into one contribution. As mentioned, this third type of interaction is reciprocal to the second type of interaction. This allows us to draw an interesting conclusion. Since an eigenmode on an infinite interconnection does not generate surface waves it is to be expected that for a semi-infinite interconnection only the part near the end, Le., near the port of the interconnection, will generate substantial surface waves.

The contribution from parts further away from the port will be less important. Each part of the semi-infinite transmission line will generate the same amount of surface waves but for the parts further away from the port there will be a more and more destructive interference between the surface waves generated by different parts. This means that we can restrict ourselves to the contribution of a small part of AI,^ and AI,^ close to A c . And from reciprocity this means that for the second type of interaction we can restrict ourselves to the contribution of distributed sources close to Ac or Bc. From the third type of interaction it is natural to deduce the fourth type of interaction. The surface waves generated by the interconnections AI,^ and AI,^ will also directly generate eigenmodes in the interconnections BI,J and B I , ~ . The circuit model for this type of interaction consists of distributed sources I B , ~ ( sd)s , I B , ~ ( sd)s , VB,~(S)ds, and VB,~(S)ds in the transmission lines which depend on the ) ,B ,I(s), and W B , ~ ( S in ) currents and voltages zA,J(S),i ~ , 2 ( s W each elementary part d s of the interconnecitions AI,^ and AI,^. Due to the fact that interconnections clo not interact very strongly with surface waves, it is clear that the first type of interaction is the most important one and the fourth type of interaction is the least important one. In this paper we will concentrate on this first type. In essence the third type of interaction is only a special case of the first type which

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

866

means that with the theory of this paper it is also possible to calculate this interaction. The second and also the fourth type of interaction can be calculated with the theory of [ 2 ] . 111. THE EXCITATION O R TRANSMITTER PROBLEM

In this section we determine the surface waves (and also space waves) generated by the component part A c when the voltages and currents V A , ~ , V A , ~ , Z A , ~and , iA,2 in the circuit model at the ports, i.e., at s = 0, are known (Fig. 1). We assume that the circuit model for the structure A was based on the reciprocity definition introduced in [2].First of all we observe that V A , J , V A , ~ , and ~ A ,Z~A ,; ~ are not fully independent, they are related through the SA-matrix. The voltage and current on both transmission lines can be written as two waves or modes, one propagating in the positive direction and one in the negative direction VA,l(s)

1

=2(Im,A,i)- exp(-jPA,is)Ki,i f 2 (I,, A,i )

iA,i(S)

exp ( j p A,i s )Ki,.j

= Im,A,i exp(-jpA,+s)Ki,;

Ac

(1)

- Im,A,i e x p ( j P A , i s ) K i , i

with i = 1, 2 and where PA,^ is the propagation coefficient of the fundamental mode on interconnection line AI;^. and K& are the excitation coefficients of the normalized fundamental mode propagating in the positive and negative direction, respectively. I m , ~ ,isi the total longitudinal current of the normalized fundamental mode. Note that the characteristic impedance Zchar,A,i is given by 2 ( I m , ~ : i ) - ’ .The model (1) is a reciprocity-current (RI) based transmission line model [ l ] for the interconnections. If E t r , ~ and , i H t r , ~ are , i the transversal field components of the fundamental mode then this mode is normalized if

$

I,/

(EtrA,i

x

H t r , A , i ) ’ 11,

dS = 1

where u, is the unit vector along the line perpendicular to the transverse plane SA,+. Using (1) it is possible to determine the excitation coefficients KA+,iand KA,i if V A , ~and i ~ ,ati s = 0 are known by simply inverting (1) for s 0. The result is 1

K A ,= ~ ;I [Im,A,iVA,i(S= 0)

- TJtVtG,(r,T’)l

. J$’(T‘>}dS‘.

(5)

The subscript ‘t’ denotes quantities parallel to the zy-plane. We restricted ourselves to the zy-components of the electric field because we will only need these components in the sequel. G A . and ~ ~ G, are scalar magnetic and electric Green’s functions, Ittis the two-dimensional (2-D) unit dyadic in the xy-plane and V t V t is a dyad. Since these Green’s functions only depend on z , z’ and the radial distance p between T and T ’ . i.e., p = 1 1 , x (T - T I ) ] , we write them as G ~ , ~ ~ ( p , z l z ’ ) and G$(T.zlz’). Due to the layered nature of the structure these Green’s functions are first determined analytically in the spectral domain and are then numerically inverse Fourier (2) transformed

S-A,%

*

where . I s . ~ . %and J t r , ~are , %the current densities of the normalized mode on interconnection AI,^ in the longitudinal and transversal direction, respectively. SA,^^ are the elements of the scattering matrix SAdefined in [2].The current density J $ ) follows from an electromagnetic simulation of the strucJf!2, and J : ) ) as in ture A. Similar current densities (Jz,!l, (4) are generated when a mode is incident along AI,^. We are now interested in the fields generated by the current densities J:), z = 1 , 2 and especially in the generated surface waves. It is possible to calculate these fields using different approaches. Since we consider planar metalizations it is suitable to use a mixed-potential formulation [3]

* 2(1m,A,z)-’iA,i(s = a)].

(3)

G ( p . z l ~ ’=) 7 ;1

l+m

G(X, zl4Jo(XP)X

(6)

where G ( p , zIz’)(G(X,zlz’)) is GA(P,zIz’)(GA(X, zlz’)) or G+(p,zJz/)(G+(X, zlz’)) and where X is the radial spectral variable. Using Cauchy’s residue theorem and after some algebra we can rewrite (6) as [4], [5]

This is the first step in the excitation problem. Given V A , ~ and i ~ ,ati s = 0 we are able to determine by what amount the eigenmodes in the lines are excited. if the structure A is excited by a normalized mode along for example AI,^ then we get the following current densities on the structure A J$,)i

= exP(jpA,ls)(-J,,A,l1I, f Jtr,A,lutr)

+ s A , 1 1 exP(-jpA,ls)(Js,A,l?~, + J t r , A , l U t r ) On (1) -

J A , ~ - S A , 1 2 exp(jpA,Zs)(Js,A,2,s

Jk“

+ Jtr,A,P%r)

AIJ

on AI,2

on A c (4)

The first term is the contribution from a branch cut, Le., the first term represents the space wave contribution. In (7) we assumed that there was a semi-infinite layer at the top of the structure and a metal plane at the bottom. if there is also a semi-infinite layer at the bottom of the structure an extra branch-cut integral has to be introduced. In (7) the wavenumber k , is given by w m with E , and pu the material parameters of the semi-infinite layer. The square roots

OLYSLACER: HIGH-FREQUENCY REClPROCITY BASED CIRCUIT MODEL

867

where the unit vector ut is in the zy-plane. The functions are the surface wave radiation patterns for the part Ac with respect to an arbitrary origin 0,i located inside A c . These functions fully characterize the far field surface wave behavior of Ac.

ft)(cp)

r *l,l

/ Po Iv. THE OBSERVATION

- - -

4

2

Fig. 4. Notations for the far field contribution of the surface waves generated by the component part A(; of structure A .

in the integrand of (7) are defined such that their imaginary part is not positive. H i 2 ) is the Hankel function of second kind and order zero with its branch-cut along the negative real axis. The second term in (7) represents the surface wave contributions. The quantities A, (Im(A,) 5 0) are the propagation coefficients of these surface waves, i.e., of the eigenmodes of the multilayered slab waveguide structure, and G,(zIz’) is the residue of G(X, zlz’) at X = X u

G,(zlz’) = lim [G‘(X,zIz’)(X X-IX”

-

Xu)].

(8)

Since the amplitude of the surface waves only decays as 1/,,@ and that of the space waves as l / r , the surface waves become dominant after some distance. This means that after some distance from Ac the main contribution in ( 5 ) comes from the surface waves. The contribution from surface wave v is given by (i) Et,u( r= ) -

2

//

I I

i

where K i , %are the yet unknown excitation coefficients of the normalized modes E t r , ~H ,% t r ,, ~in, % Ihe interconnections. For the moment we will assume that the modal part (12) of the ‘a’-field is far more dominant compared to the remaining scattered part (e.g., surface waves and space waves) at the ports of Bc. The ‘b’-field is the field generated in structure B due to an incoming normalized mode along for example the interconnection B I ,1. Along the interconnection B I , 1 the transversal field components of this ‘b’-field take the following form: Eb,tr-

AC

RECEIVER PROBLEM

V;B , I~B, , ~and , Now we want to determine the sources V B , ~ I B , at ~ the beginning of the transmission lines, i.e., at s = 0, generated by a given incident field on the component part Bc. Again a current-reciprocity transmission line model is assumed for the interconnections BI,1 and B I , and ~ the component Bc is characterized through the SB-matrix. In a first step we determine by what amount the fundamental eigenmodes in the interconnections BI,1 and Bl,? are excited through the interaction of the incident field E’”,Hi” with the component part Bc. This incident field is a field which propagates in the multilayered slab waveguide structure and consists of surface wave andlor space wave contributions. A standard approach to determine the interaction is to invoke the Lorentz recimocitv theorem. To amlv this theorem we need to define two fields, labeled and ‘b>.~ i ~ ‘a’ l , is , ~the scattered field E’‘, Hsc originating from the scattering of the incident field at the metalization of structure B.Pit the ports of Bc a part of this scattered field will consist of the eigenmodes of the interconnections B I , 1 and B1,2. The transversal components of the modal part of the ‘a’-field at these ports are given by

j ~[ ~ ~ , t t ,z ~’ () z~ ,i ~ ) ( ~ ~ p ) ~ i ) ( r(1)’ ).

OR

e x p ( j P B ,1s ) E t r , B ,1 f SB,11 expl:-jPB,

15) E t r , B ,1

+

-v,vtG+,,(z, z~)H~2)(xup).~j4)(r~)l dS’, (9) ~ t 1 L= - exp(j~B,is)~t~,B,i s D , 1 1 exp(-jPB,ls)Htr,B,l In a last step we look at the far field contribution of a surface wave. If we insert the asymptotic form of the Hankel function [6] in (9) and use the notations of Fig. 4 we obtain

(13) and along the interconnection

B1,2

E g = SB,21 e:XP(-jPB;2S)&r,B,2

~ i=:SBi ~exp ( -jpB ,2 1

2 . { GA,tt,u( z Iz’)ft’ (9) +

with fF)(p) =

X : G 4 ~ u ( z t z ’ ) [ u ‘f!)(cp)lut) t

// Ac

exp(.iX,ut

‘

ar)J2)(r’)

,2 S)jwtr,B , 2 .

(14)

In accordance with (4) the corresponding current density on Bc is denoted Jb = J g ) . Now we apply Lorentz-reciprocity (lo) theorem on the volume of which a top-view is shown on Fig. 5. This volume is bounded by the surface C which is invariant in the z-direction, i.e., C is a cylindrical surface (11) parallel to the z-axis, and which consists partly of the parts ~ the ports where of the transverse planes S B ; I , 1 and S B . I , at

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

868

B1,Z Fig. 5.

of

Cylindrical bounding surface C around the structure B c used for the application of the Lorentz-reciprocity theorem. The surface C is composed S B , ~ and , SB,Z.

Sclosingi

the modal field patterns are important and partly of a closing surface Sclosing located far enough from the structure B.In principle this closing surface can be chosen at infinity. The Lorentz reciprocity theorem takes the following form in the defined volume

//

( E , x Hb - Eb x H a ) .undS

c

=

JJ' (Eb . J ,

-

E , . J b ) dS

(15)

Bc

where un is the unit vector normal to the surface C. If we insert (12), (13) and (14) in (15) and use the fact that Eb is normal to the surface Bc we obtain

(16) It can be shown that the contribution of the closing surface negligibly small when Sclosing is located far enough from the structure B . Since the tangential component of the total electric field at the surface Bc has to be zero we have that tiz x (3'' + Ein)= 0 and we can rewrite (16) as Sclosing becomes

This expression allows us to calculate by what amount ~ excited due to the eigenmode in the interconnection B I , is interaction of the incident field with the component part Bc. It is noted that only the zy-components of E'" come in between for the considered planar structures. As mentioned earlier we had to assume that the nonmodal part of the scattered field at ~ S B , was ~ negligible. This assumption the port planes S B , and however can be dropped as is shown in Appendix B. It is shown that if the port planes are taken further down the interconnection lines the excited modes in those lines have two different origins. The first contribution is given by (17) and the second contribution comes from the interaction of the incident field with the modal currents of the interconnection lines. This last contribution corresponds with the situation of Fig. 2. In fact, Appendix B proves that it is allowed to separate the interaction with the component part from the interaction with the interconnection lines. Let us return to the situation of Fig. 1 and assume that our incident field is the field generated by Ac. In the previous derivation it is then assumed that the structure A is located outside the surface E. At the same time the scattered field has to be negligible at this surface. These two assumptions will be compatible if the effect of the presence of circuit B is negligible on the working of circuit A when circuit B is not excited at his interconnections. As a last step we have to translate the excitation coefficient K;,, to the sources V B ,and ~ I B , ~Using . the transmission line

OLYSLAGER: HIGH-FREQUENCY RECIPROCITY BASED CIRCUIT MODEL

869

AI,1

Fig. 6. Two interacting structures with arbitrarily shaped 3-D component parts and arbitrarily shaped 2-D interconnection parts embedded in a layered medium.

E2

Fig. 7. Geometry of a microstrip substrate.

equivalent defined in [2] it is easy to show that these sources are given by

V B ,1 = 2 (Im,B , 1 ) -

Kz,1

I B ,1 = Im,B,1 K z ,1 .

case, Le., (1)-(3) remain valid [l]. Let E$) and E$)denote the electric and magnetic fields in the component part of the structure A when it is excited by a normalized mode along AI,^. In fact we only need the equivalent electric and magnetic current densities .7$’ = u x Hi;) and A$) = E:) x u , on the surface Ac . These have to follow from an electromagnetic analysis of the structure A. To determine the fields generated by J $ ) and K t ) we have to generalize ( 5 ) to a coupled field mixed-potential formalism [9], [lo]

E(i)(r= )

// {

(18)

[GA(T~T’) - OOG+(rlr’)] .J$)(r’)

Ac

These two sources generate a mode in the transmission line

B I , with ~ excitation coefficient KB+,lpropagating away from Bc. If the ‘b’-field in the Lorentz reciprocity theorem is the field generated by a normalized mode incident on B c along B I , then ~ we can determine vB,2 and I B . ~ .

Ac

V. GENERALTHREE-DIMENSIONAL STRUCTURES with Up to now we concentrated on planar perfectly conducting metalization structures. In this section we will generalize G A = GA,ttItt f vtGA,taUz UzvtGA,,t f G A , ~ ~ u , ~ , the previous results to general three-dimensional (3-D), not G F = GF,ttItt f vtGF,tzUz U,vtGF,,t f G F , ~ ~ u ~ u ~ necessary perfectly conducting, objects embedded in layered (20) media. This allows the incorporation of finite conductivity and ,GA,zt, z G A ,,GF,tt, ~ ~ GF,tz,GF,,t,G F + , finite thickness which is important for MCM’s. In this way it where GA,tt,G A , ~ is in principle possible to apply the theory also to dielectric G + , and G+ are scalar Green’s functioiis for the magnetic vector potential, respectively, the electric vector potential, waveguide circuits. Consider the structure of Fig. 6 which is a 3-D gener- the electric scalar potential and the magnetic scalar potential. alization of the structure on Fig. 1. The interconnections Again all these Green’s functions only depend on p , z and AI,^, AI,^, B I , and ~ B I , are ~ now general open waveguides, z’ and hence we can proceed as in (6)-(11). However the which have been studied in for example [7], [8]. The compo- component part A c is now described by two sets of surface nent parts Ac and BC are arbitrarily shaped 3-D structures wave radiation patterns f:)( p) and g$) (p) connecting the interconnections. The same assumptions are = exp(jA,u. ~ r ) ~ $ ) d( r~’ ’) made concerning the modal character of the interconnections as in the planar case. Remark that AI,^, A I J ,A c , B I , ~B,I . ~ , Ac and Bc indicate the external surface of the structures A and B. gp) = exp(jA,u. AT) Kil)(r’)dS’. (21) First, we concentrate again on the transmitter problem. The Ac transmission line description remains the same as in the planar

+ +

ft)

// //

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

870

port B

3d

L

'

d = 3.04 mm d

P'x t A I

+3d+

dl Fig. 8. Geometry of the metalization of two square patches.

Also for the observation or receiver problem we can proceed as in Section IV. By carefully applying Lorentz reciprocity theorem as is done for the 2-D case in [2] it is possible to show that (17) has to be generalized to

Radiation pattern If (cp).u,l(Amm) 90"

Hin)is the incident field It has to be emphasised that (Ein, generated by structure A in the layered medium without the presence of structure B. VI. APPLICATION

Now we will apply the technique of Section 111 and Section IV to the dominant surface wave of a microstrip substrate. Fig. 7 shows the geometry of the substrate consisting of two layers above a perfectly conducting plate. The bottom layer has a thickness h and permittivity €1 and the top layer has infinite thickness and permittivity € 2 . Both layers are nonmagnetic, i.e., ,ul= p2 = PO. The spectral scalar Green's functions G A , ~ ~ ( and X ) G$(X) for z = z' = h are easily determined and given by

270" Fig. 9. Surface wave radiation pattern on the left patch of Fig. 8.

l u t .f 1

($)I

of the current distribution

The residue G ~ , of J G+(X) for this surface wave is

with

with z = 1, 2. The propagation coefficient A 1 of the dominant surface wave mode, which is a TM-mode, is solution of

As an example we consider the configuration of Fig. 8 at a frequency of 3 GHz. The substrate has the following characteristics: h = 3.17 ",€,,I = 11.7 and cr,2 = 1. The metalization pattern consists of two square patches with dimensions 9.12 mm x 9.12 mm, separated by a horizontal distance L H and a vertical distance L v . Both the patches are exited by a microstrip line with width d = 3.04 mm. The port of each patch is taken 3.04 mm down the microstrip lines. The characteristic impedance Zchar = 48.04 R and the propagation

OLYSLAGER: HIGH-FREQUENCY RECIPROCITY BASED CIRCUIT MODEL

coefficient p = 2.9151k0 of the fundamental mode of the microstrip lines were calculated with a classical 2-D spectral domain moment method analysis. With a commercially available electromagnetic simulator for planar structures (HPMomentum from HP EEsof) the current distribution on one patch was determined when excited at its port. The SAAparameter obtained by this simulator is given by SAA = 0.255 0.959j. At 3 GHz the propagation coefficient of the sole surface wave above cut-off is given by A 1 = 1.02226ko and the corresponding residue G ~ , J ( A = ~ ) 14075~'.Fig. 9 shows the radiation pattern I f l ( c j ) . a t / corresponding to the current distribution of the patch. Note that for the given frequency this is almost a perfect dipole radiation pattern. Fig. 10 shows the modulus of the $AB-parameter, calculated with the theory described in previous sections as a function of the horizontal distance L H when L v = 0. The dots on the figure are results obtained from an electromagnetic simulation of the total 2-port structure consisting of patch A and patch B. A good agreement is found between both results and as expected this agreement becomes even better for larger distances when the surface wave becomes more dominant. Also the phase of SAB (not shown on the figure) is in good agreement. In Fig. 11 results are presented when L v = 0.2 m. Fig. 12 shows an equivalent circuit for the whole structure of Fig. 8. The impedance 2 is given by

871

ISAB'

+

Z = 2(1+ SAA)

G ( 1 - SAA)

(27)

Fig. 10. ISasl as a function of the horizontal distance LH between the patches of Fig. 8 when LV = 0. The full line is obtained with the theory of this paper and the dots are results from an electromagnetic simulation of the global structure.

ISAJ 0.018 Surfacewave coupling Full EM-simulation

0.014

and the sources V, and I B depend on the currents i~ and are given by

~SAB v, = 1&(1iA SAA)

VA and I A follow from symmetry. VII. CONCLUSION The interaction between structures embedded in layered media was investigated. It was shown that each component can be represented by a surface wave radiation pattern. On the other hand the influence on a component of an incoming electromagnetic field was investigated by means of the Lorentz reciprocity theorem. This influence was represented by voltage and current sources at the ports of the circuit model of the component.

APPENDIXA For notational simplicity we assumed in the main text that the interconnections were monomodal. In this appendix the equations are generalized to multimodal interconnections which propagate more than one mode. Important examples of multimode interconnections are coupled microstrip lines and coupled coplanar waveguides. In the circuit model these multimodal interconnections are represented by coupled transmission lines [1].

0 0.5 1 1.5 2 2.5

3 3.5 4 4.5 5 LHb)

Fig. 11. \SAB( as a function of the horizontal distance L H between the patches of Fig. 8 when Lv = 0.2 m. The full line is obtained with the theory of this paper and the dots are results from an electromagnetic simulation of the global structure.

If the interconnection A1,i propagates 1vA,i modes then the voltages and currents on the corresponding set of coupled transmission lines are written as !!A ,i xz

iA,i

+,i (1;,A,i ) - 1 exp ( -.7g~,i')KA

+ 2@:,A,i)-

1 exP(J'p,,,is)K,,,

Lm,A,i e x p ( - j p ~ , i s ) K -t A,i

- Im,A,i e X P ( j g ~ , i s ) h ' : i , i .

(29)

This is the generalization of (1). Q ~ and , j~A , i are column matrices with N,4,i elements, & I is a column matrix with the N A , ~excitation coefficients of the modes, & A , i is a

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

872

Fig. 12. Equivalent circuit for the structure of Fig. 8.

component part A c . The remaining part of the transmitter problem is easily generalized to multimodal interconnections. We just want to mention the generalization of (11)

where 1;)is a column matrix with Ni elements representing the surface wave radiation patterns of surface wave I/ generated by each of the modes at all the interconnections. The receiver problem equations (12)-( 14) are easily generalized as in (31). Finally (17) and (18) are generalized as

and

APPENDIXB

I

(b)

Fig. 13. Close-up of the interconnection B I , to ~ demonstrate the separability between interaction with the component part and interaction with the interconnection part.

N,+ x NA,%matrix with the longitudinal currents on each line for each mode [1] and is a diagonal matrix with the propagation coefficients of the different modes. The inversion of (29) at s = 0 becomes

PA,%

K2,%= $

2(Im.A,z)-1iA,%(s = 011 .(30) If the structure A is excited by normalized modes along for example AI,^, then we get the following current densities on the structure A [ I m , A , z l A , % ( S= 0 )

In this Appendix we will show that it is allowed to separate the interaction of an externally incident field with the interconnection lines from the interaction with the component part. Let us concentrate on the close-up on Fig. 13(a) of the interconnection B I , ~ and the component part Bc.In stead of taking the port at SB,1, i.e., s = 0, we now take the port further down the line at Sb,l, i.e., s = A. In this way the component part increases by an amount ABc. The currents Jb of the ‘b’-field on ABc are modal and take the following form (see also (4)):

OLYSLAGER: HIGH-FREQUENCY REClPROCITY BASED CIRCUIT MODEL

We can rewrite this as

.lo

a

G,1= K21 +

KR+,lb)exP(jPB,ls) ds A

+SB,11i

K-&(S)

exP(-jDo.ls) ds (37)

where we used (17) for the first term on the right-hand side of (36). In the two other terms on the right-hand side we separated the longitudinal integration in the s-direction from the transverse integration. The transverse integration is incorporated in the definition of the distributed excitation coefficients K i , l ( s ) d s and K;,l(~) ds

(38) where c is a transversal integration line on ABc. The last two terms in (37) are represented in the circuit-model as distributed sources V B , ~ ( ds S ) and I B , ~ ( sds ) in the transmission line B I , ~with the theory of [ 2 ] . Remark that the longitudinal integrals of the excitation coefficients K i , l(s) ds and K B , l ( ~ )ds with the phase factors exp(jpB,ls) and exp(-jpB,ls) correspond to the coefficients P and Q defined in [2]. The first term on the other hand is represented by the discrete sources V B ,and ~ I B ,given ~ in (18). Hence, the excitation coefficient KLT, can be interpreted to consists of three parts indicated by the arrows on Fig. 13b. The first part is directly generated by the discrete source at s = 0 or in other words by the interaction of the incident field with the component part Bc. The second and third part are generated by the distributed sources or in other words by interaction of the incident field with the interconnection line. The second part is a direct wave and the third part an indirect wave which partly reflects at the component part, i.e., at s = 0.

ACKNOWLEDGMENT The author would like to thank N. FachC and his team at HP EEsof-Alphabit in Belgium and D. De Zutter for carefully reading the manuscript and for his useful suggestions.

REFERENCES [ l ] F. Olyslager, D. De Zutter, and A. T. de Hoop, “New reciprocal circuit model for lossy waveguide structures based on the orthogonality of the eigenmodes,” ZEEE Microwave Theory Tech., vol. 42, no. 12, pp. 2261-2269, Dec. 1994.

873

D. De Zutter and F. Olyslager, “High-frequency reciprocity based circuit model for the incidence of electroma.gnetic waves on general waveguide structures,” IEEE Microwave Theory Tech., vol. 43, no. 8, pp. 1826-1833, Aug. 1995. R. F. Harrington, Field Computations by Moment Methods. New York: Macmillan, 1968. L. B. Felsen and N. Marcuvitz, Radiation and Scattering of Waves. Englewood Cliffs, NJ: Prentice-Hall, 1973. S. Barkeshli and P. H. Pathak, “Radial propag;ition and steepest descent path integral representations of the planar niicrostrip dyadic Green’s function,” Radio Science, vol. 25, no. 2, pp. 161-174, Mar.-Apr. 1990. M. Abramowitz and I. A. Stegun, Eds., Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables. New York: Dover, 1970. F. Olyslager, D. De Zutter, and K. Blomme, “Rigorous analysis of the propagation characteristics of general lossless and lossy multiconductor transmission lines in multilayered media,” IEEE Trans. Microwave Theory Tech., vol. 41, no. 1, pp. 79-88, Jan. 1993. F. Olyslager and D. De Zutter, “Rigorous boundary integral equation solution for general isotropic and uniaxial anisotropic dielectric waveguides in multilayered media including, losses, gain and leakage,” IEEE Trans. Microwave Theory Tech., vol. 41, no. 8, pp. 1385-1392, Aug. 1993. J. R. Mautz and R. F. Harrington, “Electromagnetic scattering from a homogeneous material body of revolution,” Arch. Elek. Ubertragung, vol. 33, no. 4, pp. 71-80, Apr. 1979. K. Umashankar, A. Taflove, and S. M. Rao, “Ekctromagnetic scattering by arbitrary shaped three-dimensional homogeneous lossy dielectric objects,” IEEE Trans. Antennas Propagat., voll. 34, no. 6, pp. 758-766, June 1986.

Frank Olyslager (S’90-M’94) was born in Wilrijk, Belgium, on November 24, 1966. He finished high school in 1984 and received the electrical engineering degree from the IJniversity of Ghent, Belgium, in July 1989. In 1903, he received the Ph.D. degree from the Laboratory of Electromagnetism and Acoustics (LEA) of the University of Ghent with a thesis entitled Ellectromagnetic Modeling of Electric and Dielectric Waveguides in Layered Media. At present he is a Postdoctoral Researcher of the National Fund for Scientific Research of Belgium in the Department of Information Technology (the former LEA) of the University of Ghent. From 1989 until 1993 he was a Research Assistant of the National Fund for Scientific Research of Belgium. His research deals with the use of integral equation techniques to numerically solve Maxwell’s equations. His activities focus on the electromagnetic wave propagation along high frequency electrical and optical interconnections in multilayered isotropic and bianisotropic media, on the singulanty of electromagnetic fields at edges and tips, and on the study of Green’s dyadics in bianisotropic media He is also investigating the construction of transmission line models for general waveguide structures and different aspects of electromagnetic compatibility (EMC) problems on printed circuit boards and microwave circuits. He spent several scientific visits at the Electromagnetics Laboratory at the Helsinki University of Technology. He is author or co-author of about 35 papers in international journals and about 25 papers in conference proceedings He is also co-author of the book Electromagnetic and Circuit Modeling of Multiconductor Lines (Clarendon Press, 1993) in the Oxford Engineering Science Senes. In 1994, Dr. Olyslager became laureate of the Royal Academy of Sciences, Literature, and Fine Arts of Belgium In May 1995 he received The IEEE Microwave Prize. Currently he is Associate Editor of Radio Science

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

874

Higher-Order Vector Finite Elements for Tetrahedral Cells J. Scott Savage, Student Member, IEEE, and Andrew F. Peterson, Senior Member, IEEE

Abstruct- Edge-based vector finite elements are widely used for two-dimensional (2-D) and three-dimensional (3-D) electromagnetic modeling. This paper seeks to extend these low-order elements to higher orders to improve the accuracy of numerical solutions. These elements have relaxed normal-component continuity to prohibit spurious modes, and also satisfy Nedelec’s constraints to eliminate unnecessary degrees of freedom while remaining entirely local in character. Element matrix derivations are given for the first two vector finite element sets. Also, results of the application of these basis functions to cavity resonators demonstrate the superiority of the higher-order elements.

I. INTRODUCTION HE FINITE element solution of three-dimensional (3-D) electromagnetic problems using the lowest order vector finite elements, defined by Nedelec on tetrahedra [l], has been well documented [2]-[3]. These elements are commonly referred to as edge elements or Whitney elements. Because the functions do not impose normal-component continuity between cells, they produce no spurious modes in the numerical solution of the curl-curl equation. However, these mixed-order elements, which allow a constant tangential, linear normal (CTLN) representation of the fields on mesh edges, limit the accuracy of the finite element solution. Higher order basis functions, also proposed by Nedelec, allow for more accurate solutions of 3-D problems, while retaining the benefit of permitting no spurious modes. These functions fall in the general class known as “curl conforming” since they do not impose complete continuity, but do ensure tangential continuity between cells. The next higher order basis functions on tetrahedra provide a linear tangential, quadratic normal (LT/QN) representation of the fields. The basis functions of next higher order have a quadratic tangential, cubic normal (QT/CuN) representation for the fields. This article reviews these basis functions and provides closed-form expressions for the element matrices arising from the CT/LN and LT/QN functions. In addition, numerical results for the resonant frequencies of 3-D cavities are presented to illustrate the relative accuracy of the higher-order functions and the error trends as the cell sizes are reduced.

a given cell in the finite element mesh. L, is the simplex coordinate associated with node i of the cell. L, is unity at node i and decays linearly to zero at the other three nodes of the cell. For edge based functions, i and j represent the two node indices associated with that edge. For face based functions, i , j , and IC represent the three node indices at the vertices of that face. Since C T L N basis functions have six unknowns associated with any cell, they lead to 6 x 6 element matrices, while LT/QN and QT/CuN basis functions result in 20 x 20 and 45 x 45 element matrices, respectively. Also, C T L N functions have one unknown per edge throughout the global model, while LT/QN functions have two unknowns per edge and two unknowns per face. QT/CuN functions have three unknowns per edge, six unknowns per face, and three unknowns per tetrahedron. Thus, higher order basis functions lead to more unknowns for a given finite element mesh. Also, higher order basis functions result in global matrices with greater density (more nonzero entries per row and column). Therefore, the computational burden in creating and solving the finite element matrices for a given mesh increases with the order of the basis functions. The basis function definitions in Table I apply to an individual tetrahedron. Since many tetrahedra may share a certain edge, the global basis functions on that edge straddle each of those tetrahedra. Similarly, up to two tetrahedra may share a common face, so the global basis functions on that face straddle those tetrahedra. This convention ensures tangential field continuity across tetrahedra boundaries throughout the mesh.

111. ELEMENT MATRIXDERIVATIONS Efficient finite-element analysis of electromagnetic fields in 3-D regions requires computation of two element matrices associated with the curl-curl form of the vector Helmholtz operator [3]. These two matrices are

DEFINITION 11. BASIS FUNCTION Table I shows an unnormalized simplex-coordinate representation of CTLN, LT/QN, and QT/CuN basis functions for

and

Manuscript received August I , 1995; revised February 15, 1996. This work was supported in part by NSF Grant ECS-9257927 and Electromagnetic

Sciences, Inc. The authors are with the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0250 USA. Publisher Item Identifier S 0018-9480(96)03802-1.

where B, represents the ith vector basis function and V indicates integration over one tetrahedron. It is implied that (1) and (2) involve only the portion of each basis function

0018-9480/96$05.00 0 1996 IEEE

SAVAGE AND PETERSON: HIGHER-ORDER VECTOR FINITE ELEMENTS FOR TETRAHEDRAL CELLS

875

TABLE I 3-D VECTORBASISFUNCTIONS ON TETRAHEDRA (LTIQN)

(QTICuN)

6 functions

20 functions

45 functions

6 Edge Based

12 Edge Based

18 Edge Based

for all i < j ,

for all i # j ,

for all i # j ,

L,VL] - L]VL,

4VL1

& ( 2 ~- ~, ) V L , , i + j

8 Face Based

for all i < J ,

for all i < j < k , L,L1VLk-&L,VL,

&L,(VL,- V L 1 ) 24 Face Based

r,L1VLk- L]L,VLE

forall i # j # k , L,(2L, - 1)(L1VLk- LkVL,) G(L1VLk- L,VL])

which lies in a particular cell. A closed-form derivation of these matrices facilitates efficient formation of the global finite element system of equations.

TABLE I1 LOCALFACE,EDGE,AND NODECONVENTIONS Node 1

A. CT/LN Elements This section presents the derivation of element matrices for constant tangenthear normal CTLN basis functions. These basis functions are associated with tetrahedra edges and are defined in Table I as

In this representation, Lil is the simplex coordinate associated with the first node of edge i , Liz is the simplex coordinate associated with the second node of edge i, and li is the length of edge i. The simplex coordinates for a given cell are

Li = ai +biz

+ ciy +.d;z

i = l , . . .,4

(4)

and the gradient of any simplex coordinate is

V L i = b;0 + c;fj

+ d;;.

(5)

The simplex coefficients, a; . . . d; , can be computed by inverting the coordinate matrix

I

I

4

2

5

2

1 3 1

I/

4

Face ##

1

Node 1

I

Node2

4

I

2

I

3

U 1 1 Node3

I

4

I

pqi2 -

1

Fig. 1. The node and edge labeling convention used in this document. Face labeling conventions are presented in Table 11.

where (zz,yi, z i ) is the location of node i . As an intermediate step in computing the inverse of the coordinate matrix, the volume of the tetrahedron, V, is computed. As implied in (6), these simplex coefficients are not normalized to the volume of the tetrahedron.

All tetrahedra are given a local structure as illustrated in Fig. 1. The first local node associated with any edge, i l , is the lesser of the two node numbers adjacent tio the edge. Table I1 presents the convention used for local node, edge, and face numbering. For notational and computational convenience, the following notations are adopted. For two given nodes, i and j , a vector matrix and scalar matrix are defined as

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

876

~

i

=j VLi x

VLj

= ?(cidj

-

Cjdi)

+ G(bjdi

-

Using (15), (13) reduces to

b,dj )

+ q i c j - bjCi)

(7)

and 4ij

+ cicj + didj

= U L i . ULj = bibj

(8)

where each of the simplex coefficients are defined in (6). These matrices are constant for a given cell, and may be obtained as the first step in element matrix computation. Since both V and q5 are independent of position, either may be removed from any integrand. Also, note that U i j = - U j i . The evaluation of the element matrix in (1) requires the curl of each basis function

v x Bi = v x li(LilVLi2 =l i 0 x

-

Li2VLil)

(LZlVLiZ) -

1;v x

(9)

~ ~21 and~ i 2 are i the ~endpoints , of edge i , With this, (1) becomes r

The second element matrix, given in (2), requires the calculation of basis function dot products -

B; Bj -

L,2VLil)

’

lj(LjlVLj2 - LjZVLjl)

I

. v L j 2 ) - LilLjz(vLi2 . V L j i ) . VLj2) + L22Lj2(VLil VL,l) .

= lilj

This section presents the linear-tangent, quadratic-normal (LT/QN) element matrices. LT/QN basis functions exist in two forms, edge based and face based functions. The edge based functions can be written in terms of the two simplex coordinates which correspond to the endpoints of an edge, while the face based functions can be written in terms of the three simplex coordinates which correspond to the vertices of a face. The edge based LT/QN basis functions are

(Li2VLil)

= 2liVLi1 x VLi2 = 2liVil,i2.

= li(LilVLi2

B. LT/QN Elements

‘

where “el” denotes the first type of edge basis function, and “e2” denotes the second type of edge basis function. The two face-based basis functions associated with face i are

where “fl” represents the first type of face basis function, and “E” represents the second. For face-based basis functions, 21\22, and i 3 indicate the three vertex indices of face i. The element matrices for LT/QN elements involve interactions between the four types of basis functions. Therefore, the element matrices can be represented as block matrices Eelel

~ e l e 2 Eelfl

~ e l f 2

~ e 2 e l ~ e 2 e 2 ~ e 2 f l ~ e 2 f 2

(1 1)

Applying the notation of (8), (1 1) becomes

Eflel

Efle2

Eflfl

Eflf2

Ef2e1

Ef2e2

Ef2fl

Ef2f2

where, for example

Eelf1= 23

The second element matrix may then be written

Fij = l i l j [ q 5 i 2 , j 2

LiiLjl

IV

This expression may be simplified by employing the general integration formula for 3-D simplex coordinates

3!i!j!k!l! V. (3 + i + j IC I ) !

+ +

In (13), two simplex coordinates (possibly the same) are involved in each integral. These integrals can be expressed in matrix form as r2

1 1

x

R1. V x Bsl dV.

(20)

The subscript i in (20) is an edge index and the subscript j is a face index. This is implied by the superscripts, “el” and “fl”, respectively. The second element matrix, F , can be represented similarly

dV - 4i2,ji

V

.I,v

11

F=

I

Felel

~ e l e 2 pelf1

FeZel

pe2e2

~ e 2 f l pe2f2

Fflel

Ff1e2

pflfl

Ff1f2

Ff2el

Ff2e2

pf2f1

Ff2f2

]

~ e l f 2

.

(21)

To evaluate each block matrix in (19), the curl of each of the four types of basis functions is needed. For the first type of edge based basis function, “el”, the curl is

v x B,el = v x

(liLilVLi2)

= IiVLil x VLi2 = liG&iZ.

(22) Similarly, the curl of the second type of edge based basis function, “e2”, is

SAVAGE AND PETERSON: HIGHER-ORDER VECTOR FINITE ELEMENTS FOR TETRAHEDRAL CELLS

The curl of the first type of face based basis function, "fl", is slightly more complicated.

An examination of the original element matrices, (1) and (2), reveals that both are symmetric. Therefore, only those matrices on or above the main diagonal in (19) and (21) need to be evaluated. Using the curl expressions in (22)-(25) and the integration matrix notation of (15), the i , j entries in each of the submatrices in (19) are E e23l e 1 = V l i l j ( F i l , i 2 . U j l , j 2 )

(26) -Eele1

*

E:;e2 = V l i l j ( i j i l , i 2 . j J j 2 , j l ) = E e 2 e 2 = V l i l j ( V i 2 , i l . T i j 2 , j l ) = E ;ejl e l 23

(27) (28)

(30)

(31) (32)

(33)

-

. v j l , j 3+ M a. 2 . ~v.321,23 . . v3.1 .~ 2 - Mi3,jlGil,i2 * v j 2 , j 3 - 2Mi3,j2vil,i2 ' u j l , j 3 - M 2. 3 .~8321,22 . . . E3. 1 .~)2 f2f2 =V ( M i l , j l v i 2 , i 3 uj2,j3 f 2Mil,j2vi2,i3 E 2.3 . u j l , j 3 h'fZl,j3ca2,i3 . G j l , j 2

.

(34)

'

+

+

f 2MiZ,jlvil,i3 ' v j 2 , j 3 4Mi2,j2vil,i3 ' v j l , j 3 f 2Mi2,j3vil,i3 ' v j l , j 2

+-M i 3 , j l C i l , i 2 '

"

7-'j2,j3

+ 2Mi3,j2vil,i2

vjl,j3 f Mi3,j3vil,i2 ' v j l , j 2 )

(35)

877

The block entries in the second element matrix, (21), follow a similar derivation, to produce

F23f l f 2= v ( d % 3 , j 3 p i l , z 2 , j l , j 2 - & 3 , , j l p i l , i Z , j 2 , j 3 -

4. , j 3 p z. l ,.~ 22

.

3

.

~

+1'$%2,ilPil,i3,j2,.f3) ~ 2

(44)

FfZf2 = v ( h 3 , j 3 p i l , i 2 , j l , j 2 - d%3,jlPZl,iZ,j2,j3 23

- 4 .z 1 .~ p3z2,z3,31,32 . . . . + 4il,jlPi2,;3,j2,j3).

(45)

The new integration matrices, N and P., are straightforward extensions of M from (15) (46)

Iv. APPLICATION TO RESONANT CAVITY ANALYSIS The element matrices derived above may be used to construct global finite element matrices. Tlhis is accomplished by summing the element matrix terms for each tetrahedron in the mesh. With knowledge of the connectivity matrix for the mesh, it is possible to predetermine the sparsity structure of the global matrices. This allows for memory efficient construction of the global matrices in which only nonzero terms are stored. When combined with an iterative eigenvalue solution algorithm, this provides a memory and processor efficient finite element algorithm. The CT/LN and LT/QN basis functions were implemented in the finite element analysis of cavity resonators. To illustrate the relative accuracy of these functions, they were used to estimate the wavenumbers of the vector Helmholtz equation for homogenous media

The resulting matrix eigensystem was solved for the wavenumbers k , using a sparse solver based on the method of subspace iteration [4]. A rectangular cavity, with dimensions 1 x 0.5 x 0.75 m, was discretized using a commercial software package into six meshes of various density. This package creates unstructured tetrahedral meshes which strive to keep all tetrahedra well-shaped. Table I11 presents the results for CTLN basis functions and Table IV presents the results for LT/QN basis functions. The edge lengths in these tables indicate the average length of all the edges in the mesh. Although the convergence

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

878

TABLE 111 NUMERICAL RESLLTSFOR C T L N BASISFUNCTIONS 75 0.44493

0.31161

4.76396

5.10871

I

100

1

1

163

227

I

299

I 0.29486 1 0.25917 1 0.23571 1 0.21848 I Wavenumber 5.17303

5.14063

5.19660

5.19896

5.23599

TABLE IV NUMERICAL RESULTS FOR LT/QN BASISFUNCTIONS

of any particular eigenvalue is usually erratic, the general error behavior of the two methods is apparent when the average error of several modes is visualized. Fig. 2 shows the average percent error of the first eight modes plotted versus the average length, h, of all edges in the mesh. A curve fit through the data points indicates the order of convergence of the two methods. For this geometry, a curve fit to the C T L N data has an exponent of 1.98, while a similar curve through LT/QN data has exponent 3.86. This O ( h 2 )convergence for C T L N elements and O(h4) convergence for LT/QN elements is consistent with previous 2-D numerical investigations and the theoretical 3-D dispersion analysis of Warren [ 5 ] . Similar trends were observed for other cavity geometries including spherical and cylindrical shapes. As indicated in Fig. 2, on a given mesh, the LT/QN basis functions give much more accurate solutions for the wavenumbers, IC, than do C T L N basis functions. However, for a given mesh, LT/QN basis functions require more unknowns than C T L N functions. Based on experience, approximately six times as many unknowns are needed for the LT/QN basis functions on a given mesh. Even though LT/QN elements require more unknowns, Fig. 3 demonstrates that for equal numbers of unknowns, LT/QN elements still outperform C T L N elements. Higher-order elements also affect the global matrix sparsity characteristics. C T L N basis functions typically yield global matrices with approximately 15 nonzero entries per row, while LT/QN functions usually give global matrices with approximately 35 nonzero entries per row. This means that LT/QN solutions require more computational effort than

3-D Error Comparison Average of First 8 Modes 10

................................

-

I

0.01

0.2

0.3

0.4

0.5

Average Edge Length, h Fig. 2. Error comparison between CTLN and LT/QN elements based on the average mesh edge length, h. The slope (exponent) of the CTLN curve fit is 1.98, while the slope of the LT/QN curve fit is 3.86.

CT/LN solutions with an equal number of unknowns. However, in analyses with equal unknowns, C T L N basis functions require a mesh with many more tetrahedra. This increased mesh density increases the mesh creation time and the mesh storage requirements. Therefore, LT/QN basis functions used with relatively coarse meshes provide more efficient solutions. This assumes that no additional geometry modeling errors are introduced by using a coarse mesh. For arbitrary geometries,

879

SAVAGE AND PETERSON: HIGHER-ORDER VECTOR FINITE ELEMENTS FOR TETRAHEDRAL CELLS

REFERENCES

3-D Error Comparison Average of First 8 Modes

J. C. Nedelec, “Mixed finite elements in R3,” Num. Math., vol. 35, pp.

315-341, 1980. J. F. Lee and R. Mittra, “A note on the application of edge-elements for modeling three-dimensional inhomogeneously-filled cavities,” IEEE Trans. Microwave Theory Tech., vol. 40, pp. 1161-1113, Sept. 1992. J. Jin, “The finite element method in electroimagnetics.” New York: Wiley, 1993 Y. Li, S. Zhu, and F. A. Fernandez, “The efficient solution of large sparse nonsymmetric and complex eigensystenis by subspace iteration,” ZEEE Trans. Magn., vol. 30, pp. 3582-3585, Sept. 1994. G. S. Warren, “The analysis of numerical dispeirsion in the finite-element method using nodal and tangential-vector elements,” Ph.D. Dissertation, Georgia Institute of Technology, Atlanta, Georgia, 1995.

10

10

Inn

innn

inono

Unknowns Fig. 3. Error comparison between CT/LN and LT/QN elements based on the number of unknowns. The curves diverge, indicating that LT/QN elements provide more accurate solutions with equivalent computational effort.

this necessitates the use of curved tetrahedra, a subject of future work.

V. CONCLUSION This paper has reviewed the first three mixed order, 3-D, vector basis functions for tetrahedra. Element matrices were derived for C T L N and LT/QN elements. Results were presented from the application of these two basis function sets to cavity resonator analysis. It is concluded that higher order basis functions are preferable to lower order basis functions, since higher order bases are capable of providing more accurate results with coarse tetrahedral meshes, fewer unknowns, and less overall computation. Numerical solutions for the resonant wavenumbers were observed to converge at an O ( h 2 ) rate for the C T L N functions and an O(h4) rate for the LT/QN functions.

J. Scott Savage (S’95) received the B.E.E. degree from the Georgia Institute of Technology in June, 1992, and the M.S.E.E. degree from the University of Kentucky in December, 1993. Since 1994, he has been a Ph.D. student and ]Research Assistant at the Georgia Institute of Technology. From 1992 to 1993, he ‘was a Research Assistant studying crosstalk in cable bundles with Ford Motor Company. His interests include computational and applied electromagnetics. Mr. Savage is a member of Eta Kappa Nu and Tau Beta Pi.

Andrew F. Peterson (S’U,-M’83-SM’92) received the B.S., M.S., and Ph.13. degrees in electrical engineering from the University of Illinois, UrbanaChampaign in 1982, 1983, and 1986, respectively. From 1987 to 1989, he served as Visiting Assistant Professor at the University of Illinois. Since 1989, he has been a member of the faculty of the School of Electrical and Computer Engineering at the Georgia Institute of Technology, where he is now an Associate Professor. He teaches electromagnetic field theory and computational electromagnetics, and conducts research in the development of computational techniques for electromagnetic scattering, microwave devices, and electronic packaging applications. He is an Associate Editor of the IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION. Dr. Peterson is the immediate past Chairman of the Atlanta chapter of the IEEE AP-SMTT society. He is also a Director of ACES and a member of URSI Commission B, ASEE, and AAUP.

880

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

Riccati Matrix Differential E Formulation for the Analysis of Nonuniform ultiple Coupled Microstrip Lines Jen-Tsai Kuo, Member, IEEE

the received signal. If the lines are electrically short and made on a low-loss substrate, the radiation can dominate the loss mechanism. Reflections are caused by the position-dependent impedance values along the lines. Note that all the aforesaid factors depend on the nonuniformity of the lines and on the operation frequency which make the characterization of the NCML’s network become a complicated task. Several methods have been developed to analyze multiple NCML’s. Mehalic and Mittra [4] investigated the tapered multiple microstrip lines using a spatial iteration-perturbation approach technique. Oh and Schutt-Aine [5] analyzed the nonuniform lines based on a time-domain scattering parameter formulation incorporated with the closed-form expressions of voltage variables for divided short uniform lossless lines. Mao and Li proposed a method of convolution-characteristic [6] and I. INTRODUCTION a method of equivalent cascaded network chain [7] to handle ONUNIFORM coupled microstrip lines (NCML’s) play the transient response of NCML’s. Palusinski and Lee [3] an important role in both analog and digital microwave used Chebyshev polynomials to expand the current and voltage integrated circuits. Using NCML’s, for example, a folded allalong the nonuniform lines in the time-domain to predict the pass two-port network [ l ] and a directional coupler [l], [2] reflections and cross talk of general multiple coupled line can be realized with high coupling values operating over systems. an ultra-wide frequency band. To date, NCML’s serve as In frequency domain, Arabi et al. [8] presented an electrithe interconnections in most chip packages for digital intecal field integral equation formulation based on a combined grated circuits of switching speed covering the microwave or approach of using closed-form near and far field approximamillimeter-wave regime [3]-[9]. With the advances of today’s tions for the Sommerfeld microstrip Green’s functions. The semiconductor fabrication technology, the major portion of accuracy of this technique can be set to any desired value. delay time in a microwave integrated circuit (MIC) can be In [9],Pan and his colleagues extended the method in [3] to due to these interconnection lines [3]. One possible way for the frequency domain. The advantages of analyzing NCML’s reducing the delay time is to increase the density of the in the frequency domain over the time domain were also interconnecting NCML’s. As the NCML’s become shorter discussed. or are placed closer, the nonuniformity of the lines must be To calculate the input reflection coefficient matrix of terproperly designed in order to obtain transmitted signals with minated NCML’s, we derive a differential matrix equation, sufficiently high quality. which is known as the Riccati matrix differential equation When high-speed signal travels along NCML’s, the received (RMDE), based on a reciprocity-related definition of the line signal at the load end can be degraded due to 1) dispersion, 2) voltages and currents for hybrid-mode multiple coupled micross talk, 3) losses, and 4) reflections. The cross talk and discrostrips. The RMDE is expressed in terms of the normal mode persion are due to the differences of relative effective dielectric parameters of coupled microstrips and solved by method of constants for different modes and at different frequencies, moments. The method of solution is also extended to calculate respectively. The losses which include conductor, dielectric, the scattering parameters of 2N-port NCML’s networks. and radiation attenuation factors will lower the power level of The presentation is organized as follows. Section I1 describes the background of the mathematical modeling of Manuscript received August 22, 1995; revised February 15, 1996. This work was supported in part by the National Science Council, Taiwan, Grant NCML’ s and lists the mathematical formulas to describe the NSC 85-2213-E-009-002. reflection along the lines. Section I11 presents the method of The author is with the Department of Communication Engineering, National solution to the nonlinear RMDE. In Section IV, the converChiao Tung University, 1001 Ta-Hsueh Road, Hsinchu, Taiwan, R.O.C. gence behavior of the analysis method is investigated and Publisher Item Identifier S 0018-9480(96)03783-0.

Abstruct- A Riccati matrix differential equation (RMDE) is formulated for analyzing nonuniform coupled microstrip lines (NCML’s) in the frequency domain. The formulation is based on a reciprocity-related definition in the theory of multiconductor transmission lines under quasi-TEM assumption. The hybridmode nature of modal phase velocities and strip characteristic impedances for multiconductor microstrip structure is included. The nonlinear RMDE is first transformed into a first-orderlinear differential matrix equation which can be efficiently solved using method of moments. A convergence study is performed to investigate the sufficient number of basis functions used in the method. The voltage-scattering parameters of a tapered microstrip and two three-line structures are presented. The frequency responses of a pair of nonuniform coupled lines are measured and compared with calculated results.

0018-9480/96$05.00 0 1996 IEEE

~

KUO: RICCATI MATRIX DIFFERENTIAL EQUATION FORMULATION

88 1

Fig. 1. A system of N-conductor nonuniform coupled microstrip lines.

several numerical aspects are discussed. Numerical results for certain nonuniform single microstrip and. three-line structures are presented and discussed. Section V compares the measured frequency responses of a nonuniform two-line structure with the calculated results. Finally, Section VI draws the conclusion. 11. THERICCATIMATRIXDIFFERENTIAL EQUATION (RMDE)

It is known that a system of uniform N-conductor coupled microstrip lines and a ground line support N dominant or quasi-TEM modes. For the NCML’s in Fig. I, we neglect the fringing fields, which produces radiation loss, caused by the gradual change of waveguide cross section. For lines with abrupt discontinuities, field-theoretical oriented formulations, such as that in [SI, can be referred to enhance accuracy of results. At any z along the NCML’s, through the full-wave solution, an N x N matrix [MI], called the eigencurrent matrix, can be obtained [lo]. Of [MI] each column vector consists of total currents on the lines for a given mode. Based on the orthogonality of modal voltage and current vectors, an eigenvoltage matrix [Mv] iis uniquely defined [ I l l . For each mode, inner product of the eigenvoltage and eigencun-ent vectors is set to be the total electromagnetic power transfer. This is an important fact that lleads our field problem to be able to be formulated by circuit quantities. The characteristic admittance matrix along the lines is given by [111

The equivalent distributed capacitance matrix [C] and inductance matrix [L]along the lines can be derived [lo]

[CI = [MI]diag ( P k I W ) [Mv] [LI = [Mvl diag ( P k I W ) [&‘I]-l

where w is the angular frequency and p k iis the phase constant of the kth mode. Note that all the entries in [MI], [Mv], [Yc], [L], and [C] are dispersive and position-varying along the NCML’s. Let [I]and [VI be the line current and line voltage column vectors of which the kth entries are the total current and voltage on the kth line, respectively. Then, from the multiconductor transmission-line theory [3]

[I]= [Y,nI[Vl [VI’ = -PI [I1 [I]’ = -[YI[VI

(4) (5) (6)

where [Kn]is the input admittance matrix :seen at z toward the load, [Z] and [Y] are, respectively, the series impedance and shunt admittance matrices per unit length of the NCML’s, and the prime (’) represents the derivative with respect to z . If the tapered lines are lossless, [Z] = j w [ L ] and [Y] = jw[C]. Let the reflection coefficient matrix along the longitudinal direction be [ p v ] , then [Y,,] and [pv] are related by [lo]

[Knl = [YCl([Ul - [Pvl)([Ul + [pvl)r1

(7)

where [U] is the identity matrix of size N x N . Substitution of (4) and (5) into (6) leads to

[LI -’ [Knl[~l[Y,,l+ [ Y ]= 0. It can be shown that [Yc] is symmetric and the important aspect of reciprocity is guaranteed. If the load network has an admittance matrix identical to the [Yc] of the NCML’s at the load end, then there is no reflection.

(2) (3)

(8)

Inserting (7) into (S), one obtains [PVI’ = j([YI[Pvl + [PVI[‘Yl)

+ ([UI + [Pvl)[c:I([ul - b v l )

(9)

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

882

Ma where [r]= [Yc]-'[Y] = [Z][Yc] = [ M ~ ] d i a g ( , & ) [ - M \ - ] - ~ 1x1= C[XmIc,,,(z) (16) . that (8) and (9) are known as and [GI = [ Y ~ ] - ~ [ Y c l ' / 2Note m=O the RMDE [12] which is nonlinear. It is believed that [13] is the first literature that formulated the RMDE (8) for studying where C, (2) is the shifted Chebyshev polynomial of order m of the first kind defined over 0 5 z 5 L, L being the length general nonuniform transmission lines. In the case of single nonuniform line, (9) becomes the of the NCML's, and [X,] and [A,] are constant matrices. Riccati scalar differential equation (RSDE) [ 121. To simplify The matrix [A] is first expanded into a linear combination this nonlinear differential equation, many authors [ 14-16] of C m ( z ) of which each coefficient matrix can be obtained neglected the p i term. The solution of pv at z = 0 can by Gaussian Chebyshev quadratures [19]. Then [Am]s can be then be obtained through a simple Fourier transform of a obtained since each C m ( z )is known as a polynomial of z of function of the line characteristic impedance. Based on the degree m. Following the method in [20],more precisely the transform, synthesis of matching transformers and couplers Galerkin procedure in method of moments, one can obtain using nonuniform transmission lines have been developed [ 151, [XI [Ql[W)I (17) [16]. Note that the legitimacy of the negligence is relied on the fact that p: and Fig. 9 shows the results of E:. These results are compared with the theoretical complex permittivity values as given by the Cole-Cole equation

ing integral equation derived by Levine and Papas [15]:

& r*

where dielectric constant of probe filler; intrinsic admittance of the probe; radius of the outer conductor: b U radius of the inner conductor; E complex permittivity of sample; Ydut measured input admittance of the probe when brought in contact with sample under test; J , (x) cylindrical Bessel function of argument 2 ; and where

=&l

r -JET

!I

E,

Yo

2T k=-. (3) A, If the values for the standard material are not found in literature then the standard material is measured using cavity perturbation techniques (over limited frequency range) and those values are assumed to be the theoretical values. The corrected admittance of the material under test Y d u t c is then found by using Ydvtc

= (7corrYdut

(4)

where C,,,, is the correction factor obtained using measurements on the standard material, and Y d u t is the measured admittance of the material under test. Once the corrected value of the admittance is found then the complex permittivity results ( E : ) can be obtained by substituting the values of into (2). For our purpose, we used a Mathcad@ subroutine to carry out these calculations.

where ionic conductivity of the liquids; relaxation time; a distribution parameter; E, optical relative permittivity; E, static relative permittivity; W radian frequency. The Cole-Cole parameters have been taken from the National Bureau of Standards data [15]. For the methanol measurements the calibration of the probe was done by taking a response calibration with respect to short and using deionized water as the standard material. Figs. 10 and 11 show the obtained E’, results of ceramics that have been measured at room temperature with the metallized ceramic probe. For Fig. 10 the material under test was A1203 + 0.1 mol% MgO, and the standard material used for the measurements was ZrOz 8 mol% Y2O3. For Fig. 11 the material under test was ZrOz 8 mol% Y203. and the standard material was A103 0.1 mol% MgO. The results in both cases are compared to the previously mentioned commercially available HP probe, and to cavity perturbation results [16]. When a line is fit to the data points it shows very good correlation between the three methods. However there are a few erroneous points in the metallized probe results. This is due to the more relaxed calibration method. It may be U

T

+

+

+

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

930

25

20 13 I

Theoretical Probe

23 12

4

i

21 19

x i

15

0.5

1.0

1.5

2.0

2.5

1

1

0

o

1 4

>

0.0

3.0

0.5

1.0

Fig. 9. E: of methanol measured with metallized probe versus theoretical values for the frequency range from 500 MHz to 3 GHz.

2.0

2.5

3.0

+

Fig. 11. E: of ZrOz 8 mol% Y 2 0 3 measured after a response calibration versus cavity techniques and the commercially available Hewlett-Packard probe.

I

j

m Met. Probe H P Probe o Cavlty

0 cavity I Water

21 E,’

9

19

’1

5

1

0.5

1.0

1.5

2.0

2.5

”

(Stndrd)

o Alumina (Stndrd)

23

0.0

1.5

Frequency (GHz)

Frequency (GHz)

15

Met. Probe HP Probe Cavity

Ii 1

3.0

Frequency (GHz)

0.5

1.0

1.5

2.0

2.5

3.0

Frequency (GHz) Fig. 10. E:! of A1203 +O. 1 mol% MgO measured after a response calibration versus cavity techniques and the commercially available Hewlett-Packard probe.

Fig. 12. E: of ZrOz 8 mol% Y 2 0 3 measured after a full calibration using two different standard materials versus cavity techniques.

shown that in a calibration where only two standards are used that the results will not be as accurate as a method that uses the more complete calibration procedure. Figs. 12 and 13 show the results of a more complete calibration on ZrOz 8 mol% YzO3. For these measurements the full 12-term error correction model was employed by using the three standards needed for the calibration. An open, a short, and then two different standard materials were used as the standard material so that a comparison could be made. A1203 + 0.1 mol% MgO was used as one of the standard materials and deionized water was used as the other. Fig. 12 shows the results for E ; , and Fig. 13 shows the results for E:. In both figures a comparison is made between cavity perturbation techniques and measurements done on the same

material, but just merely changing the standard material used in calibration. It is clear that these results are more accurate than the results shown in Fig. 11 because of the improved calibration procedure used in this case. It can also be pointed out that the standard material does not have to be in the general dielectric range of the material under test when a more complete calibration in performed. After verification that the metallized probe worked well at room temperature the probe was used to make hightemperature measurements. Due to oxidation problems of the metallization at high temperatures, the probe must be used in a controlled atmosphere. Therefore a heating system has been developed that utilizes a molybdenum coil to heat the sample and the probe in a hydrogen atmosphere. The coil is typically

+

+

~

BRINGHURST AND ISKANDER OPEN-ENDED METALLIZED CERAMIC COAXIAL PROBE

93 1

27

-I 1

0.8

0.6 E;’

OW4 0.2

--Cavity +Water (Stndrd) - c t A l U m l n a (Stndrd)

25

I

0

23

19

-

17

0.0

0.5

1.0

1.5 2.0 2.5 Frequency (GHz)

0

B

23

Probe Cavity Probe Cavity

21 E,’

I

1.5 2.0 2.5 Frequency ( G H r )

1.0

1.0

1.5

2.0

2.5

3.0

Frequency ( G H r )

3.0

i

3l

IProbe (800 C)

Cavlty Probe 0 Cavity 0 Probe n caviiy A Probe A CaVlty 0

Er)) 2

19

I

Fig. 15. E’, of Zr02 + 8 mol% Y 2 0 3 versus cavity results at both 600 and 80OOC in the frequency range from 500 MHz to 3 GHz.

(400 C )

(400 C) (200 C) (200 C)

I

0.5

3.0

Fig. 13. E: of ZrOz + 8 mol% YzO:) measured after a full calibration using two different standard materials versus cavity techniques.

0.5

cavity ( 8 0 0 C)

0.5

(800 (600 (600 (400 (400

C) C) C) C) C) (200 C) (206 c )

1.5 2.0 2.5 Frequency (GIHz)

1.0

3.0

Fig. 14. E’, of ZrOz + 8 mol% Y z 0 3 versus cavity results at both 200 and 4 O O O C in the frequency range from 500 MHz to 3 GHz.

Fig. 16. E: of ZrO2 $ 8 mol% Y z 0 3 versus cavity results at 200,400, 600, and 80OOC in the frequency range of 500 MHz to 3 GHz.

4-in long, with a 2-in diametlx. The coaxial probe is placed inside the coil with the coil centered at the interface between the probe and the sample. Figs. 14-16 show the results of Zr02 8 mol% Y203. For these measurements A1203 0.1 mol% MgO was used as the standard material. These measurements were done with the simpler response calibration. Fig. 14 shows the results of E’, compared to cavity perturbation results at 200 and 400°C. Fig. 15 shows the results of E’, compared to cavity perturbation results at 600 and 800OC. Fig. 16 shows the results of E: compared to cavity perturbation results at 200, 400,600, and 800°C. A line has been lit to the data points in all three figures and it is clear that there is good agreement between the two measurements procedures. It is also clear that there are a few points that are not as accurate and, once again, this may be

due to the simpler calibration method used. A full calibration would undoubtfully give more accurate overall results.

+

+

IV. ERRORANALYSIS It has been mentioned earlier that a small air gap between the probe and the sample under test can cause significant errors in the measurements [5]. An analysis was therefore made to quantify the errors caused by the differential thermal expansion which is inherent in metal probes. In the heating up cycle the outer conductor will heat up faster than the inner conductor, and will therefore expand imore than the inner conductor, thereby causing a measurement air gap between the inner conductor and the material under test. In the cooling down cycle the outer conductor will cooll faster than the inner

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

932

TABLE 1 COMPARISON BETWEEN THE FDTD REFLECTION COEFFICIENT RESULTS VERSUS AKALYTICAL DATA[ 5 ] FOR THE CASE OF AIR GAPSBETWEEN THE PROBE ARE U = 3.0 MM. b = 7.0MM, AND THE PROBE FILLER HAS &: = 2.54, AND THE AND THE SAMPLE UNDER TEST. PROBE DIMENSIONS MATERIAL UNDER TESTHAS THE RELATIVE COMPLEX %RMITTIVITY E : = 10.0 - j1.0.RESULTS WERE CALCULATED AT f = 1 GHz

.o

F b e liftoff

Ir

I r (Phase in Degrees)

(Magnitude)

FDTD Model

Reference [5]

I

Reference [5]

FDTD Model

1 0.0mm

10.950

10.952

1-32.0

1-29.0

10.5 mm

10.991

10.993

1-10.8

I -10.4

I 1.0 mm

10.998

10.996

I -7.1

1-7.0

I I I

TABLE I1 COMPARISON OF MAGNITUDE AND PHASE OF INPUTIMPEDAKCE WHEYAIR GAPSDUETO DIFFERENTIAL THERMAL EXPANSIONS ARE PRESENT Air gap between inner material under 0 25 mm

0 25 mm

mm

mm

mm

mm

Maqnitude

TABLE I11 COMPARISON OF COMPLEX %R'vIITTIVITY RESULTS W H E K CHARACTERISTIC IMPEDANCE VALUES ARE CHANGED DUE TO DIFFERENTIAL TRANSVERSE THERMAL EXPANSIOK BETWEEN THE INNER AND OUTERCONDUCTORS OF THE PROBE

% difference in:

Change in characteristic impedance by:

,

+5%

-5% I

-10%

+I 0%

1

I

E,'

I

1.67%

I

2.48%

I

4.5%

I

2.76%

I

I

E.II

I

45%

I

123%

1

73.9%

I

157%

I

conductor, and will therefore shrink faster, thereby causing an air gap between the outer conductor and the material under test. Errors created by these air gaps were modeled using finitedifference time-domain (FDTD) simulations. The probe was first modeled flush against the material under test, and the input impedance was measured. Then an air gap, of 1.0, 0.5, and 0.25 mm, between the inner conductor and the material under test was modeled, and an air gap, of the same sizes, between the outer conductor and the material under test was also modeled. In all cases the input impedance was measured and compared to the original case. To check the accuracy of our FDTD model we first compared the FDTD results with those published earlier [5] for the case when the entire probe (both inner and outer conductors) were separated from the sample by an air gap. Table I shows comparison between the FDTD results and the analytical data published earlier [5] where it may be seen that there is a good agreement. Table I1 gives the FDTD results of the error analysis for the high temperature case. Table I1 shows a comparison of both the magnitude and phase of the calculated input impedance when the probe is in contact with the sample versus when it is assumed that there is a differential thermal expansion between the inner and the outer conductors as described earlier. It should be noted that the error analysis presented in this paper are different from those available in literature [ 5 ] .

Unlike available analysis where uniform gaps between the probe and the sample are examined [SI, the presented results show possible errors due to differential thermal expansions which are highly relevant to the topic of this paper. It may be desirable to report error analysis results for air gaps with dimensions smaller than 0.25 mm. This was however, very difficult using our presently available uniform mesh FDTD code. Memory requirements on the available work stations was simply prohibitive. It is possible to obtain such results using the variable mesh code presently under development [ 181 and these results will be reported in a future article. Another analysis was also completed to examine the effect of thermal expansion on the characteristic impedance of the probe and hence its effect on the accuracy of the dielectric properties measurements. The ratio of the inner and outer conductors was changed by 3 ~ 5 %and *lo%. *5% change in the characteristic impedance corresponds to k change in the radius of the outer conductor b or & change in the radius of the inner conductor a. The new characteristic impedance ( Z o )of the probe was then calculated. The admittance (Yo) can then be found, and substituted into (2) to find the new complex permittivity values. The new complex permittivity values were then compared to the original values obtained for a probe with an characteristic impedance of 50 0. Table I11 gives the results of the comparisons.

BRINGHURST AND ISKANDER: OPEN-ENDED METALLIZED CERAMIC COAXIAL PROBE

933

9

++Mag. and Phase +Phase -0- Magnil ude

8

7

+Magnltude

-- 6

2 W'

200

.-E

E

'L 5

b L 4

0 L

44

150 100

2 50 1

0

2

1

3

4

5

1

Fig. 17. Error analysis of

E',

by changing Sll of the material under test.

From the results in Tables I1 and I11 it may be noted that longitudinal and transverse differential thermal expansion of the probe may have significant effect on the dielectric properties measurements. The development of the metallized ceramic probe, therefore, is highly important and certainly needed. Finally an analysis was made to estimate how the errors in the measurement of ,911 (since the ,911 term is what is measured) could affect the complex permittivity results. An actual measurement was made and then 5'11 was intentionally varied to see how the results for the complex permittivity would change. The magnitude and phase of ,911 were changed from 1-5% (both separately and then together) and then substituted back into the 12-term error correction model to The value of the complex permittivity is obtain a new Ydut. then obtained by using (2). Figs. 17 and 18 show the absolute value of the resulting errors in the value of the complex permittivity by measurement errors in ,911 when the probe was in contact with the material under test. Similar analysis was done by altering ,911 from 1-5% (both magnitude and phase separately and then together) when calibration was made on the standard material. Figs. 19 and 20 show the errors in E: and of the material under test as a result of errors in S l l when the probe was in contact with the standard material. By examining Figs. 17-20 it can be seen that for only a 5% error in the measurement of S11, of the material under test, the error can be as high as 8.7% for E:, and as high as 289.1% for e:. It can also be seen that for only a 5% error in the measurement of ,911, of the standard material, the error can be as high as 7.1% for e:, and as high as 261.5% for E:. This leads to the conclusion that for accurate dielectric measurements using the open-ended coaxial probe, errors in ,911 For standard materials, and for the material under test should be less than 1%. Other errors that may be introduced from differential thermal expansions are expected to further reduce the accuracy of the results.

3

2

5

4

% Change in Sl,

% Change in S,l

Fig. 18. Error analysis of :E by changing S11 of ithe material under test.

+Mag. and Phase +Magnitude -0- Phase 6

.-Ew-

5

L

E 4 L u1

$ 3 2 1

0

2

1

3

4

5

% Change in S,, Fig. 19. Error analysis of

E',

by changing

S11

of the standard material.

It is shown that the use of the metallized ceramic probe will significantly enhance the accuracy and capability of this dielectric properties measurement system.

v.

SUMMARY AND CONCLUSION

An open-ended metallized ceramic coaxial probe has been developed to provide accurate dielectric properties results over a broad frequency band and up to 1000°C. The probe was made of alumina and metallized with 3-mil thickness of molymanganese and a 0.5-mil protective coating of nickel. This metallized ceramic probe ensures more accuracy than any type of metal probe because of the minimal thermal expansion difference between the inner and outer conductors. Typical frequency range of the developed probe is from 500 MHz to 5 GHz. At lower frequencies, a larger coaxial probe is expected to provide improved accuracy, while a coaxial probe with

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

934

+Magnltude

=w-

/ I

200

.-C

2

1

3

4

5

% Change in SI, Fig. 20. Error analysis of

E!

by changing

SI1 of the standard material.

smaller inner and outer conductor dimensions, are expected to help improve the accuracy of these measurements at higher frequencies. Two different calibration methods have been employed in the presented measurement results. One method consists of a response calibration with either a short-circuit or an opencircuit, and a standard material, which should have similar complex permittivity values as the material under test. The other calibration method is the three calibration standard method used to determine the unknown coefficients in the 12term error correction model in the vector network analyzer. There is a trade-off in deciding which calibration to use. The former method is much simpler and less time consuming than the later, but the results are not as accurate. In the later method any standard material can be used as the calibration standard as long as the dielectric properties are well known. An error analysis has been given to show that relatively small errors in the measurement of 5’11 of either the material under test or the standard material can result in significant errors in the complex permittivity results. The resulting errors are much greater in E: than in the E; values. This is due to open-ended coaxial probes inherent problem of not being suitable for measuring the loss factor accurately in particularly low-loss materials. The broadband capabilities of this probe, however, and the simplicity of use continues to make this measurement technique an attractive and frequently used one. An error analysis has also been given to show that the differential thermal expansion which is inherent in metal probes can lead to erroneous results in the measurement of the complex permittivity. This can be due to both air gaps between the probe and the material under test, and changes in the characteristic impedance of the inner and outer conductors of the probe. REFERENCES [I] E. C. Burdette, F. L. Cain, and J. Seals, “In vivo probe measurement technique for determining dielectric properties at VHF through mi-

crowave frequencies,” IEEE Trans. Micro. Theory Tech., vol. MTT-28, pp. 4 1 4 4 2 7 , 1980. E. Tanabe and W. T. Joines, “A nondestructive method for measuring the complex permittivity of dielectric materials at microwave frequencies using an open transmission line resonator,” IEEE Trans. Instrum. Meas., vol. IM-25, pp. 222-226, 1976. D. Misra er al., “Noninvasive electrical characterization of materials at microwave frequencies using an open-ended coaxial line: Test of an improved calibration technique,” IEEE Trans. Micro. Theory Tech., vol. 38, pp. 8-14, 1990. 0. M. Andrade, M. F. Iskander, and S. Bringhurst, “High temperature broadband dielectric properties measurement techniques,” Microwave Processing of Materials Ill, vol. MRS-269, pp. 527-539, 1992. J. Baker-Jmis, M. D. Janezic, P. D. Domich, and R. G. Geyer, “Analysis of an open-ended coaxial probe with lift-off for nondestructive testing,” IEEE Trans. Instrum. Meas., vol. 43, pp. 711-718, 1994. (Errata for the dimensions of the probe listed in the captions of Figs. 4 and 5 as well as for mislabeling Figs. 10 and 11 will be published by the authors of [5] soon). P. De Langhe, K. Blomme, L. Martens, and D. De Zutter, “Measurement of low-permittivity materials based on a spectral-domain analysis for the open-ended coaxial probe,” IEEE Trans. Instrum. Meas., vol. 42, pp. 879-886, 1993. C. Li and K. Chen, “Determination of electromagnetic properties of materials using flanged open-ended coaxial probe-full-wave analysis,” IEEE Trans. Instrum. Meas., vol. IM-44, pp. 19-27, 1995. D. Blackham, “Calibration method for open-ended coaxial probe/vector network analyzer system,” Microwave Processing of Materials III, vol. MRS-269, pp. 595-599, 1992. M. F. Iskander, and J. B. Dubow, “Time-and frequency-domain techniques for measuring the dielectric properties of rocks: A review,” J. Microw. Power, vol. 18, no. 1, pp. 55-74, 1983. R. M. Hutcheon et ai., “RF and microwave dielectric measurements to 1400’C and dielectric loss mechanisms,” Microwave Processing of Materials III, vol. MRS-269, pp. 541-551, 1992. S . Bringhurst, “Metallized ceramic open-ended coaxial probe for broadband high temperature dielectric properties measurements,” M.S. thesis, Univ. Utah, Salt Lake City, UT, Mar. 1995. S. Bringhurst, M. F. Iskander, and 0. A. Andrade, “New metallized ceramic coaxial probe for high-temperature dielectric properties measurements,” Microwaves: Theory and Application in Materials Processing II, Ceramic Trans., vol. 36, pp. 503-510, 1993. S. Bringhurst, M. F. Iskander, and P. Gartside, “FDTD simulation of an open-ended metallized ceramic probe for broadband high-temperature dielectric properties measurements,” Microwave Processing of Materials IV, vol. MRS-347, pp. 221-228, 1994. I. Fitzpatrick, “Error models for systems measurements,” Microwave J., vol. 21, no. 5, pp. 63-66, 1978. H. Levine and C. H. Papas, “Theory of the circular diffraction antenna,” J . Appl. Phys., vol. 22, no. 1, 1951. Tables of Dielectric Dispersion Data f o r Pure Liquids and Dilute Solutions, National Bureau of Standards Circular 589, Nov. 1, 1958. S. Bringhurst, M. F. Iskander, and 0. M. Andrade, “High-temperature dielectric properties measurements of ceramics,” Microwave Processing of Materials I l l , vol. MRS-269, pp. 561-568, 1992. Z. Huang, J. Tucker, and M. F. Iskander, “FDTD modeling of realistic microwave sintering experiments,” in IEEE AP-S Symp. Dig., vol. 3, pp. 1798-1801, 1994.

Shane Bringhurst (S’95) received the B.S. and M.S. degrees in electrical engineering from the University of Utah, Salt Lake City, UT, in 1992 and 1995, respectively. He is currently working on a Ph.D. degree at the University of Utah in the area of electromagnetics and microwaves. He is both a Research and Teaching Assistant. His research includes: dielectric properties measurements, microwave sintering of ceramics, microwave and RF drying of ceramics, and numerical techniques in electromagnetics.

BRINGHURST AND ISKANDER: OPEN-ENDED METALLIZED CERAMIC COAXIAL PROBE

Magdy F. Iskander (F’93) received the B.Sc (Hon.) degree from the University of Alexandria in 1969, and M.Sc. and Ph.D. degrees from the University of Manitoba, ‘Winnipeg,Man., Canada, in 1972 and 1975, respectively. He is Professor of Electrical Engineering at the University of Utah, Salt Lake City. He is also the Director of the NSFIIEEE Center for Computer Applications in Engineering Education (CAEME) and Director of the State Center of Excellence for Multimedia Education and Technology. In 1986, he established the Engineering Clinic Program to attract industrial support for projects to be performed by engineering students at the University of Utah. Since then, more than 60 projects have been sponsored by 21 corporations from across the United States. He is also the Director of the Conceptual Learning of Science I:COLOS) USA Consortium which is sponsored by Hewlett-Packard Company and has eleven Member universities from across the United States. He edir.ed a book on Microwave Processing of Materiuls (Materials Research Society, 1994), and authored a textbook on Electromagnetic Fields and Waver, (Englewood Cliffs, NJ: PrenticeHall, 1992). He edited the CAEME !ioftwure Book, vol. I, 1991; Vol. 11, 1994; and co-edited two other books on Microwave Processing of Materials, (Materials Research Society, 1991 and 1992). He edited two special issues of the Journal of Microwave Power, trne on “Electromagnetics and Energy Applications,” March 1983, and the other on “Electromagnetic Techniques in Medical Diagnosis and Imaging,” September 1983. He also edited a special issue of the ACES Journal on computer-aided electromagnetics education and the Proceedings of the 1995 International Conference on Simulation in Engineering Educution. He is the editor of the journal Computer Applications in Engineering Education (CAE), (New York: Wiley). He has published more than 150 papers in technical joumals and made numerous presentations in technical conferences. His present fields of interest include the use of numerical techniques in electromagnetics. Dr. Iskander has received the Curtis W. McGraw ASEE National Research Award for outstanding early achievements, the ASEE George Westinghonse National Award for innovation in Engineering Education, and the 1992 Richard R. Stoddard Award from the IEEE EMC Society. He is a Distinguished Lecturer for the Antennas and Propagation Society of IEEE. He is a member of the National Research Council Committee on Microwave Processing of Materials.

935

936

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

DC Instability of the Series Connection of Tunneling Diodes Olga Boric-Lubecke, Member, IEEE, Dee-Son Pan, Member, IEEE, and Tatsuo Itoh, Fellow, IEEE

Abstruct- An oscillator with a series connection of tunneling diodes produces significantly higher power than a single diode oscillator. However, a circuit with series-connected tunneling diodes biased simultaneously in the negative differential resistance (NDR) region of the I-V curve is dc unstable. This dc instability makes the series connection oscillator fundamentally different from a single diode oscillator. Associated with the dc instability are the phenomena of minimum oscillation amplitude and frequency. Due to the minimum oscillation amplitude, it is critical to provide the impedance match hetween the oscillator circuit and the series connection at the desired oscillation amplitude level. An in depth, comprehensive analysis of the dc instability is given here. Based on this analysis, a numerical procedure is developed to accurately predict the minimum oscillation amplitude and frequency. Time domain simulations which give further insight into series-connection oscillator behavior are discussed. The effect of increasing the number of diodes on the oscillator performance is explored as well. Based on numerical and simulation results, oscillators with several tunnel diodes connected in series were designed and tested. Experimental results that confirm the existence of the minimum oscillation amplitude are presented for oscillators with two, three, and four tunnel diodes.

I. INTRODUCTION resonant tunneling diode (RTD), currently the fastest room temperature solid-state active device, is considered to be a promising millimeter- and submillimeter-wave source. However, RTD oscillators reported so far have not produced useful power levels. The maximum power generated by an RTD oscillator at microwave frequencies to date is 20 mW at 2 GHz [I]. At submillimeter frequencies, 0.2 pW at 420 GHz with a GaAs/AIAs diode [2], and 0.3 pW at 712 GHz with an lnAs/AISb diode were reported [3]. Only if the power levels generated by these diodes are increased, will RTD’s be useful in practical applications. Besides the fundamental thermal and impedance constraints [4j, the output power from a single RTD oscillator is also limited by stability considerations [5]. To meet typical system requirements, it would be necessary to combine the output power from several RTD’s. Several power-combining schemes have been proposed for oscillators using tunneling diodes (the term “tunneling diode” will be used to refer to an RTD or a p-n tunnel diode). For example, a modification of the Kurokawa-Magalhaes combiner was used Manuscript received October 17, 1995; revised February 15, 1996. This work was supported by the Joint Services Electronics Program Grants AFOSR F49620-92-C-0055 and F49620-96-C-0002. 0. Boric-Lubecke is with NASA/JPL, Pasadena, CA 91 109 USA. D.-S. Pan and T. Itoh are with the Department of Electrical Engineering, University of California at Los Angeles, Los Angeles, CA 90024 USA. Publisher Item Identifier S 0018-9480(96)03796-9.

to combine the power from two RTD oscillators at 75 GHz [6j. The parallel connection of 25 RTD’s was successfully used to generate 5 mW at 1.18 GHz [7]. Quasi-optical power combining was also proposed based on an RTD oscillator with a slot-coupled quasi-optical open resonator [8]. A 16-element tunnel diode grid oscillator successfully operated at 2 GHz [9]. The series connection of tunnel diodes in order to increase the oscillator output power was proposed and successfully demonstrated at low frequencies in 1965 by Vorontsov and Polyakov [lo]. The series integration of RTD’s in order to enhance the output power of an RTD oscillator at millimeterwave frequencies was proposed by Yang and Pan [ 111. An oscillator with several tunneling diodes connected in series produces significantly higher power than a single diode oscillator, but that is not the only difference between the two oscillators [IO], [ll]. Due to the negative differential resistance (NDR) region in the dc I-V curve of a single tunneling diode, a circuit using several tunneling diodes biased simultaneously in the NDR region and connected in series is dc unstable. This means that if there is no RF signal present in the circuit, the diodes cannot stay biased simultaneously in the NDR region. However, if there is an RF signal present, the bias points may be maintained simultaneously in the NDR region (the circuit may be RF stable), provided that the RF signal satisfies certain conditions. Owing to the dc instability, the design of an oscillator with a series connection is much more involved than just determining the impedance of several diodes in series. A simple dc battery is insufficient to bias several tunneling diodes simultaneously in the NDR region. Associated with the dc instability are the phenomena of minimum oscillation amplitude and frequency. Due to the minimum oscillation amplitude, it is critical to provide the impedance match between the oscillator circuit and the series connection at the desired oscillation amplitude level. An in depth, comprehensive analysis of the dc instability is given in this paper for the piece-wise continuous, linear approximation of the tunneling diode I-V curve. The existence of minimum amplitude and frequency was recognized in [lo], and their values were estimated numerically in [11]. In this paper, the minimum oscillation amplitude and frequency are physically explained and analytically derived based on the dc instability analysis. A numerical procedure is developed to accurately predict these two parameters. MWSpice simulations which give further insight into series-connected tunneling diode oscillator behavior are discussed. The effect of increasing the number of diodes on oscillator performance is explored as well. Based on numerical and simulation results, oscillators

0162-8828/93$03.00 0 1993 IEEE

~

BORIC-LUBECKE et al.: DC INSTABILITY OF THE SERIES CONNECTION OF TUNNELING DIODES

931

with several tunnel diodes connected in series were designed and tested. Experimental results that confirm the existence of minimum oscillation amplitude are presented for oscillators with two, three, and four tunnel diodes.

11. DC INSTABILITY In the following analysis, we will consider two tunneling diodes connected in series and biased with one dc battery [Fig. l(a)]. First we assume an RF signal present in the circuit, and analyze the stability of the difference between the total voltages on individual diodes (diode operating points). Next we analyze the dc bias distribution if there is no RF signal present. Each diode is represented as a parallel connection of a capacitance and a voltage controlled current source. A series resistance is not taken into account to simplify the analysis. Diodes are assumed to have the same capacitance and identical I-V curves. For the purpose of this discussion, a piece-wise continuous linear approximaticin of the dc I-V curve may be used [Fig. l(b)]. In the NDR region, the slope of the I-V curve is -R;', where R i l :is positive. For simplicity, we can assume that the slope of the I-V curve in the positive differential resistance (PDR) regions, R;', is the same on the first and second rising branches. Ideally, since the diodes are identical, the total applied dc and RF voltage should be equally divided between the two diodes. However, due to device noise [ 121, there is always a small difference between individual diode voltages, AVd(t)

avd(t)

Vdl(t)- VdZ(t).

(1)

/

-4 (b)

Initial amount of noise AVd(0) can be estimated from the device physics avd(0)

a&,.

(2)

If there is no RF signal present in the circuit, the total instantaneous voltage will be equal to the applied dc bias voltage

+

Vdl(t) Vd2(t)= Vdc.

+

= o.

Vdl(t) and Vd2(t)can be obtained, which using the initial condition (5) give solutions

(3)

Differentiating (3), we obtain dVdl(t)

Fig. 1. Simplified model of an oscillator with two tunneling diodes in series (a) and a piece-wise continuous linear approximation of the de I-V curve of a single tunneling diode (b).

(4)

dt dt In this case, we can assume that the initial bias distribution is

For both diode operating points on the first or on the second rising branch (R,) or in the NDR region (R,), the differential equation for the total current reduces to (Fig. 1)

Substituting (1) into (6), we obtain the differential equation for Avd, which using the initial condition (2) gives the solution

(7) If there is no RF signal present, (3) and (4) can be substituted into (6). Rearranging the terms, differential equations for

Therefore, if both diode operating points are on the first or on the second rising branch, any difference between diode voltages due to noise will decrease proportionally to the diode RpC constant (7). If there is no RF signal, the bias voltage in the range Vdc 5 2V, and Vdc 2 2V, will be equally divided between the diodes. In this case the bias distribution is dc stable, since any initial difference decreases in time and eventually vanishes (8). Any time two or more diodes are biased simultaneously in the NDR region, if no RF signal is present but there is a difference in bias voltages due to noise, this difference will grow proportionally to the diode R,C constant (8). Given that both diodes were initially biased close to the middle of the NDR region, the diode bias points will switch to the PDR region when the difference in bias voltages becomes equal to the extent of the NDR region. Therefore, having both diodes biased simultaneously in the NDR region is an unstable solution for the bias distribution. However, if an RF signal is present, bias points can be maintained simultaneously in the NDR region.

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

938

If one diode’s operating point is on the first rising branch, and the other on the second rising branch, and if the diode bias points are in the NDR region, the diodes must be oscillating out of phase. However, it is impossible for the diodes to oscillate out of phase and have the same total current. Therefore, in this case, the difference between diode operating points will become the difference between dc bias points. Assuming that the operating point of diode one is on the first rising branch and of diode two on the second rising branch, differential equation for the total current becomes

C -dvdi dt

vdl

f -

R,

dVd2 xC-+-

dt

vd2

Rp

-Ib.

(9)

60.-

-

-

4

E

-

I

i

-

-

Voltage [VI

(a)

Using (3) and (4), and rearranging the terms, differential equations for Vdl(t) and V&(t) can be obtained. Assuming that diode one is initially biased at the peak voltage, and diode two at the valley voltage, we can find solutions given by

Since the quantity in brackets in (10) is always positive [Fig. l(b)], V& ( t ) will decrease, while Vd2(t) will increase until they reach the equilibrium values. At this point conduction currents are equal. If the equilibrium voltages found from (10) are perturbed, solving (9) gives the solutions

Therefore, this bias voltage distribution is dc stable, since it returns to the initial value if perturbed. However, if an RF signal is applied externally to the diodes, and if the slopes of the first and second branch of the dc I-V curve are nonlinear and not equal, this bias distribution may be RF unstable [13], ~41. This analysis may be extended for additional diodes connected in series. Furthermore, assuming a nonlinear dc I-V curve, implies that R, and Rp are not constant, and the differential equations must be solved numerically. However, the same qualitative conclusions as for the piece-wise continuous dc I-V curve case will still hold. The biasing problem can be explained from this dc stability analysis. As the bias voltage is increased slowly from zero to 2Vp,both diodes will be biased on the first rising branch, and bias voltage will be equally divided between the diodes. As soon as diodes are simultaneously biased in the NDR region, if there is any difference in individual bias voltages AVd due to noise, this difference will start growing. Effectively, bias voltage will be increasing on one diode, and decreasing on the other diode (7). Hence, the rate of the increase of the bias voltage must be greater than the rate of increase of AVd, so that the voltage at each tunnel diode may be increased as well. Otherwise, if the rate of increase of bias voltage is slow compared to the diode R,C time, the diode bias points will switch to the PDR region. If a dc bias voltage sufficient to bias all tunneling diodes in the middle of the NDR region is applied gradually, the dc instability will divide this voltage

Voltage [VI (b)

Fig. 2. DC I-V curve of (a) three vertically integrated RTD’s for 50 pm diameter device and (b) two series connected tunnel diodes (back diodes MIX1 168 manufactured by Metelics Co.), measured using an HP 4145 curve tracer. Due to the dc instability, multiple current peaks are observed.

so that all the diodes are biased in the PDR region. The dc I-V curve of the series connection exhibits multiple peaks, because the diodes cannot simultaneously be biased in the NDR region. Fig. 2 shows the dc I-V curve of (a) three vertically integrated RTD’s for 50 pm diameter device and (b) two series connected tunnel diodes (back diodes M l X l l 6 8 manufactured by Metelics Co.), measured using an HP 4145 curve tracer. There are several effective solutions to the biasing problem: fast electric pulse excitation [lo], RF excitation 1131, [ 141, optical illumination [15], and successive triggering [lo], ~41. Besides the biasing problem, oscillation amplitude and frequency limitations are important consequences of the dc instability as well. Both of these phenomena stem from keeping diode bias points simultaneously in the NDR region during oscillation, as will be explained in the following sections. 111. MINIMUMOSCILLATION AMPLITUDE A signal generated by a single tunneling diode oscillator may be of any amplitude smaller than the maximum oscillation amplitude (V,f),,, at which the negative differential conductance becomes zero [ 131. For a series connection oscillator, if the oscillation amplitude is so small that the diode operating points remain in the NDR region during the whole oscillation period, the difference between diode voltages will be steadily increasing (7). When this difference becomes comparable to

BORIC-LUBECKE et al.: DC INSTABILITY OF THE SERIES CONNECTION OF TUNNELlNG DIODES

the extent of the NDR region, diode bias points will switch to the PDR region and the oscillatijon will cease. If the oscillation amplitude is large enough so that during each oscillation period diode operating points cross to the PDR region, oscillation may be maintained. During the time diode operating points are in the PDR region, the difference between diode voltages decreases (7) canceling some of the increase that occurred in the NDR region. To maintain the oscillation, the oscillation amplitude must be large enough to cover a portion of the PDR region sufficient to cancel the growth of the difference between diode voltages, AV,, during each period (7)

939

diode operating points are in the NDR region, the difference in diode voltages AVd increases (7). If this difference becomes comparable to the extent of the NDR region, the diode bias point will switch to the PDR region and the oscillation will cease. During one oscillation period, the NDR region is swept twice, and the time spent in the NDR region continuously without crossing to the PDR region is one half of the total time spent in the NDR region, t,. Therefore, using (7), the maximum time that can be spent in the NDR region during one oscillation period, can be found from the following equation

AVd o e ( t n L x 1 2 R n C = v -. where tn and t, are the times spent in the NDR and PDR region respectively during one oscillation period. For better accuracy, nonlinear I-V curves can be divided into small segments where a linear I-V relationship can be established. Assuming that -G,,, is the slope of ith segment in the NDR region, and G,, the sum of the slopes of j t h segments in the two PDR regions, (12) can be simplified as k

m

1

1

where k is the known number of segments in the NDR region, and rrL is the number of segments in the PDR regions to be determined from (13). When (13) is satisfied, minimum oscillation amplitude can be calculated as

where AV is the voltage segment. Equation (14) was derived for the steady-state oscillation ( 12), without taking into account the amount of noise that accumulates during oscillation buildup. In reality, minimum amplitude will be somewhat higher, and it will slightly depend on 1 he oscillator configuration and type of excitation. For relatively low frequencies, minimum oscillation amplitude is also related to the oscillation frequency, as will be shown in the following section. For a single tunneling diode oscillator, due to a broad range of values of negative resistance and the absence of a low frequency limitation, an oscillation is likely to occur even if impedance matching is not very accurate, but output power may be very low. However, in the case of an oscillator with several tunneling diodes in series, without appropriate impedance matching, oscillation is not possible at all. Due to the minimum oscillation amplitude, it is critical to provide the impedance match between the oscillator circuit and the series connection at the desired oscillation amplitude level. IV. MINIMUMOSCILLATION FREQUENCY For a single tunneling diode oscillator, oscillation is possible at any frequency below the high frequency cutoff fc (frequency at which the real part of the diode impedance becomes zero [4], [5]). In the case of a series-connection oscillator, there is also a lower limit on the oscillation frequency. At lower frequencies, oscillation period is longer, and diode operating points spend more time in the NDR region. During the time

vP .

(15)

For nonlinear I-V, divided into small linear segments, (15) can be modified as

1

where (&)ma,

= kA(tn)rnax.

(17)

To find the minimum oscillation frequency, we also need to know the ratio of the time spent in the NDR region to one oscillation period. Oscillation amplitude determines this ratio, and for each amplitude a corresponding minimum frequency can be found. As the oscillation amplitude increases, more time is spent in the PDR region and a smaller portion of the period is spent in the NDR region, so the minimum frequency will be lower. The maximum oscillation amplitude (VT,)max determines the lowest of the minimum frequencies. For (V,,),,,, the minimum ratio of time spent in the NDR region to one oscillation period can be found (Fig. 3)

Using (16)-( 18), the minimum oscillation frequency can be determined as

fmm

=

(tn)max

.

(19)

On the other hand, at very low frequencies, the oscillation frequency will determine the minimum oscillation amplitude. If the oscillation period is comparable to (tn)max, (12) is not valid any more, and (15) should be considered instead. With oscillation frequency known, and with (tn)maxdetermined from (17), the minimum oscillation amplitude can be determined using the same reasoning as for (18) (Fig. 3). Even though multiple NDR regions exist in the dc I-V curve of the series connection, there is a minimum oscillation frequency for the series-connected tunneling diodes. This is because the NDR region, that exists effectively during oscillation, is not what is seen in the dc I-V curve. If all diodes are oscillating simultaneously, the I-V curve of a single diode is effectively stretched N times along the voltage axis,

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

940

0.6

7

0.4

E

v

c

e

2 u

0.2

0.0

-/

0.0

1.0

1.5

2.0

2.5

Voltage [VI

0.0

VP

0.5

Fig. 4.

RTD dc I-V curve and ninth-order polynomial fit.

TABLE I MINIMUM OSCILLATION AMPLITUDE AND FREQUENCY FOR SERIES CONNECTION OF TUNNEL DIODESAND RTD’S

< Fig. 3. Oscillation amplitude determines the ratio of the time spent in the NDR region to one oscillation period.

and there is only one NDR region, N times wider than the NDR region of a single diode. Such a “stretched” I-V curve could be observed only if the voltage of a curve tracer was swept at the frequency higher than the minimum oscillation frequency. During oscillation, a small portion of a rectified, stretched I-V curve can be measured [14].

V. SIMULATION RESULTS Minimum oscillation amplitude and frequency were calculated for one tunnel diode and one RTD based on the dc instability analysis (14), (19). A low peak current tunnel diode (back diode), MIX1168 manufactured by Metelics Co., is described in [13] (diode’s dc I-V in Figs. 2 and 3), and its I-V curve was modeled using the fifth-order polynomial fit. The RTD considered here is an actual diode fabricated at UCLA, described in [16] (diode’s dc I-V in Figs. 2 and 4), and its I-V curve was modeled using the ninth-order polynomial fit (Fig. 4). All calculations were done using the Math Works Inc. MATLAB program [17], and results for both diodes are presented in Table I. Initial amount of noise was estimated as Shot noise [12], for bias current in the middle of the NDR region and a bandwidth of 200 MHz, which gave 0.04 mV for the tunnel diode and 0.03 mV for the RTD. The series-connected diode oscillator behavior was simulated in the time domain using HP-EEsof‘s Microwave SPICE program [18]. Simulations were done for oscillators with two,

three and four tunnel diodes and for oscillators with three RTD’s connected in series. The large signal impedance of each diode was calculated using the procedure described in [16], and planar microstrip oscillators were designed at 2 GHz using HP-EEsof‘s Touchstone program [ 191. Circuit impedance was matched to the device impedance at the oscillation amplitude levels above the minimum value determined from the MATLAB calculations. Oscillators were designed as two port networks, with one port used for the excitation and the other for the output. It was initially verified that the oscillation occurred at the design frequency and amplitude level. Then, circuit parameters were varied to decrease the oscillation amplitude to the value below which oscillation could not be maintained. Fig. 5(a) shows the voltage on each tunnel diode in a two-diode oscillator for the oscillation amplitude below the minimum value. Diodes are initially biased in the NDR region and are oscillating, but after about forty periods the bias points switch to the PDR region [one to the first rising branch (V&) and the other to the second rising branch ( V d 2 ) of the dc I-V curve], and oscillation ceases. Fig. 5(b) shows the voltage on each tunnel diode during stable oscillation (these voltages are equal), with the oscillation amplitude above the minimum value. As the number of diodes was increased in the same oscillator configuration, it was found that the minimum oscillation amplitude increased somewhat, from 0.125 V for a two-diode oscillator to 0.13 V for a fourdiode oscillator. Minimum oscillation amplitudes determined from the MWSpice simulation for oscillators with a series

~

BORIC-LUBECKE et nl.: DC INSTABILITY OF THE SERIES CONNECTION OF TUNNELING DIODES

94 I

0.6 1

I

0.3

E

p

-L

0.2

-5

M 0

8

c

0.1

0.0

02

nn 0 Y.”

0

10

20

30

-10

40

60

80

lo0

Fig. 6. Excitation voltage and voltage on each tunnel diode for an RF frequency of 50 MHz, and an RF amplitude of 0.125 V on each diode, for a two-diode oscillator.

(a1

0

20

Time [ns]

Time [nsl

0.0

0.4

B

20

30

Time [ns]

(b)

Fig. 5. Voltage on each tunnel diode for the oscillation amplitude: (a) below the minimum value and (b) above the minimum value for a two-diode oscillator.

connection of two tunnel diodes and three RTD’s are shown in Table I. These values agree with values obtained from the MATLAB calculations within 5%. Similar simulations were done to determine the minimum oscillation frequency for the series connection. Diodes were initially biased with a slow pulse, in the PDR regions, and an external RF signal was applied ~tothe oscillator. The frequency of the external RF signal was lowered until diode voltage was not sinusoidal any more. In the case of a two tunnel diode oscillator, for frequencies between 2 GHz and 80 MHz, and an RF voltage amplitude of 0.125 V on each diode, diode bias points were not seen to switch to the PDR region. As the frequency was lowered, the voltage amplitude had to be increased to keep diodes biased in the NDR region. Fig. 6 shows the excitation voltage (dashed line) and voltage on each tunnel diode (solid lines, in a two-diode oscillator for an RF frequency of 50 MHz, and an RF amplitude of 0.125 V on each diode. During the first period of the external RF signal, diode bias points switched to the NDR region, and diode voltages started to follow the sinusoidal external voltage. However, before the end of the second period the diode bias points switched to the PDR region. This was repeated during each consecutive period, with each diode sometimes switching to the low voltage and sometimes to the high voltage. For a 20 MHz RF frequency, the required voltage amplitude on each diode was already 0.2 V, ,which was above the maximum amplitude for which the diode can be used in an oscillator [ 131. Therefore, free-running oscillations at frequencies below

20 MHz cannot exist. In the frequency range between 20 MHz and 80 MHz, oscillations can exist, but only with an oscillation amplitude above the very large minimum amplitude, which is strongly dependent on frequency. Spurious oscillations at frequencies below 20 MHz may exist simultaneously with the desired oscillation signal of much higher frequency, but their frequency and amplitude would also be constrained by the dc instability, which makes them less likely to occur. Using the same procedure, it was determined that the minimum oscillation frequency for a series connection of RTD’s is about 50 MHz. Minimum oscillation frequency values obtained from the MWSpice simulations are about three times lower than that obtained from MATLAB calculations (Table I), due to the different initial amount of noise, which is critical for this estimate (15). In MWSpice simulations, only thermal noise in a very narrow bandwidth is taken into account, whereas a more realistic prediction of Shot noise is assumed in MATLAB calculations.

VI. EXPERIMENTAL RESULTS Due to the high series resistance and the unsuitable diode configuration, it was not possible to use the RTD’s in the experiments. Previously described tunnel diodes were used for proof-of-principle experiments. Several two-, three-, and four-diode oscillators were designed at 2 GHz for oscillation amplitudes above the minimum values determined from the simulations. The phase shift between the devices due to the diode package was taken into account for the impedance calculations [20]. Oscillators were fabricated in a one-port configuration, and RF excitation was used to trigger the oscillation [ 131, [ 141. The performance of two-diode oscillators was described in detail in [14]. Neither oscillation nor switching of bias points was observed in the oscillators designed for oscillation amplitudes smaller than 0.14 V. Therefore, it was determined that 0.14 V was the minimum oscillation amplitude, which is about 10% higher than that predicted by the simulations. This discrepancy is reasonable, taking into account that experiment must be “noisier” than the simulation. Also, in the simulation, the diodes were assumed to be identical, whereas it was impossible to avoid some variation in their dc I-V characteristics in the experiment.

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

942

-E

5

-20 -40

QI

B -60

2

-80

1.915

2.015

2.115

Frequency [GHz] (a)

0.00

0.05

0.10

0.15

0.20

Oscillation Amplitude [VI Fig 8. Possible frequency and amplitude ranges of oscillation for a series -303,dBm

I

I

-40

-60

I

I

I

I

I

I

I

I

1

I

I

1

connection of tunnel diodes.

I

k.OdGHz{

-20

I

I

I

I

I

I

1

-80

1.915

2.015

2.115

Frequency [GHz] (b)

Fig. 7. Output spectrum of the oscillator with four tunnel diodes for an excitation frequency of (a) 1.99 GHz and (b) 2.05 GHz.

Oscillators with three and four diodes were much more difficult to trigger than two-diode oscillators [20]. For oscillation amplitudes below 0.154 V neither switching nor oscillation was observed. Switching of bias points was occurring, but stable oscillation was not possible for the amplitude of 0.154 V in the case of three- and four-diode oscillators. For the oscillation amplitude of 0.176 V, the signal for the three- and four-diode oscillators existed only for a very limited bias voltage range and the output power was barely increased as compared to the two-diode oscillators. For the best three-diode oscillator, the oscillation frequency was 2.001 GHz, with an output power of -27.5 dBm, whereas a comparable two-diode oscillator gave -28 dBm at 1.995 GHz. Some increase in the minimum oscillation amplitude with the increased number of diodes was predicted from the simulation, but not such a drastic difference in the oscillator performance. This could be attributed to the diode separation due to the package, which makes the impedance of the series connection difficult to predict. Also, when the length of the series connection becomes significant compared to the circuit dimensions, generated power is no longer combined strictly on the device level. The four-diode oscillator designed for an oscillation amplitude of 0.176 V exhibited multimode operation characteristic for circuit level power combining [21]. In this case, the oscillator output signal was dependent on the excitation frequency. For an excitation frequency of 1.99 GHz, excitation was possible with power as low as -37 dBm. After the excitation signal was turned off,

the oscillation frequency was 1.976 GHz, with an output power of -29.5 dBm [Fig. 7(a)]. For an excitation frequency of 2.05 GHz, much larger excitation power (-7 dBm) was required. In this case, the oscillation frequency was 2.015 GHz, with an output power of -27 dBm [Fig. 7(b)]. It is possible that only three diodes oscillated in a four diode oscillator, especially in a lower frequency mode. A series integrated device, such as proposed in [11], should not be affected by these problems, since the diodes can be placed very close to each other. Therefore, connecting more than two packaged diodes in series does not provide further insight into the behavior of a series integrated device. It was not possible to experimentally determine the minimum oscillation frequency in the way it was done in simulations, because it would be very difficult to determine the voltage amplitude on each diode during excitation, In the simulation, we can easily see at which frequency the voltage amplitude on each diode exceeds (VTf)max.In the experiment, it is very difficult to accurately determine how much power is delivered to the diodes during excitation, and the corresponding voltage on each diode, because diode impedances are changing during excitation. Possible frequency and amplitude ranges of oscillation for a series connection of tunnel diodes are presented in Fig. 8. For a single diode oscillator, an oscillation is possible with any amplitude between 0 V and 0.182 V, provided that the frequency of oscillation is below the high frequency cutoff for that oscillation amplitude. For a series connection, stable oscillation is possible only for the oscillation amplitude above the minimum value, and in a limited frequency band (shaded area in Fig. 8). This makes the oscillator design more challenging, since large signal impedance of the series connection must be determined accurately for successful oscillator operation. VII. CONCLUSION The dc instability of the series connection of tunneling diodes and its consequences, such as the minimum oscillation amplitude and frequency, were analyzed in detail. Based on the dc instability analysis, a numerical procedure was developed to estimate the minimum oscillation amplitude and

BORIC-LUBECKE ef al.: DC INSTABILITY OF THE SERIES CONNECTION OF TUNNELING DIODES

frequency, and these parameters were calculated for one tunnel diode and one RTD. MWSpice simulations were carried out to give further insight into series-connected diode oscillator behavior. The minimum oscillation amplitude obtained from MATLAB calculations agrees within 5% with results obtained from MWSpice simulations for both diodes. Based on these results, oscillators with two, three, and four diodes were designed and tested. As observed in simulations, the minimum oscillation amplitude increased with the number of diodes in the experiment. Also, three- and four-diode oscillator performance was significantly different from the two-diode oscillator performance, possibly because of the phase shift between diodes due to the package. The four-diode oscillator exhibited a multimode operation characteristic for the circuit level power combining. Therefore, connecting more than two packaged diodes in series does not provide further insight into the behavior of a series integrated device, such as proposed in [ll]. The minimum oscillation amplitude makes the seriesconnection oscillator design more challenging than a single diode oscillator design. For a single tunneling diode oscillator, due to a broad range of valuer; of negative resistance and the absence of a low frequency limitation, an oscillation is likely to occur even if impedance matching is not very accurate, but output power may be very low. However, in the case of an oscillator with several tunneling diodes in series, without appropriate impedance matching, oscillation is not possible at all. Due to the minimum oscillation amplitude, it is critical to provide the impedance match between the oscillator circuit and the series connection at the desired oscillation amplitude level.

ACKNOW.LEDGMENT The authors would like to thank C. W. Pobanz for useful suggestions regarding the manuscript.

REFERENCES S. Javalay, V. Reddy, K. Gullapali, and D. Neikirk, “High efficiency microwave diode oscillators,” Elec. Lett., vol. 28, no. 18, pp. 1699-1700, Aug. 1992. E. R. Brown ef al., “Oscillations8up to 420 GHz in GaAs/AIAs resonant tunneling diodes,” Appl. Phys. Lett., vol. 55, no. 17, pp. 1777-1779, Oct. 1989. E. R. Brown et al., “Oscillation:; up to 712 GHz in InAs/AISb resonanttunneling-diodes,” App. Phys. Lett., vol. 58, pp. 2291-2293, May 1991. S. M. Sze, Physics of Semiconductor Devices, 2nd ed. New York: Wiley, 1981. C. Kidner, I. Mehdi, J. R. East, and G. I. Haddad, “Power and stability limitations of resonant tunneling diodes,” IEEE Trans. Microwave Theory and Tech., vol. 38, pp. 864-872, July 1990. D. P. Steenson et al., “Demonstration of power combining at W-band from GaAs/AlAs resonant tunneling diodes,” in Proc. 5th Int. Symp. Space THz Tech., Ann Arbor, MI, May 1994, pp. 756-767. K. D. Stephan et al., “5 mW parallel-connected resonant-tunneling diode oscillator,” Elec. Lett., vol. 28, no. 15, pp. 1411-1412, July 1992. K. D. Stephan et al., “Resonant-tunneling diode oscillator using a slotcoupled quasioptical open resonator,” Elec. Lett., vol. 27, no. 8, pp. 647-649, Apr. 1991. M. P. DeLisio et al., “A 16-element tunnel diode grid oscillator,” in Proc. 1995 IEEE AP-S Int. Symp., Newport Beach, CA, June 1995, pp. 1284-1287.

943

Y. I. Vorontsov and I. V. Polyakov, “Study of oscillatory processes in circuits with several series-connected tunnel diodes,” Radio Eng. Electron. Phys., vol. 10, pp. 758-763, May 1965. C. C. Yang and D. S. Pan, “Theoretical investigations of a proposed series integration of resonant tunneling diodes for millimeter-wave power generation,” IEEE Trans. Microwave Theory and Tech., vol. 40, pp. 434441, Mar. 1992. J. 0. Scanlan, Analysis and Synthesis of Tunnel Diode Circuits. Aberdeen, U.K.: WileyKJniversity Press, 1966. 0. Boric-Lubecke, D. S. Pan, and T. Itoh, “RF excitation of an oscillator with several tunneling devices in series,” IEEE Microwave and Guided Wave Lett., vol. 4, no. 11, pp. 364-366, Nov. 1994. __, “Fundamental and subharmonic excitation for an oscillator with several tunneling diodes in series,” IEEE Trans. Microwave Theory and Tech., vol. 43, pp. 969-976, Apr. 1995. 0. Boric-Lubecke and T. Itoh, “Optical illumination of series integrated resonant tunneling diode,” in Proc. 1993 IEEE AP-S and URSI Radio Science Meet., Ann Arbor, MI, June/July 1993, p. 96. 0. Boric-Lubecke, D. S. Pan, and T. Itoh, “Large signal quantumwell oscillator design,” in Proc. 23rd Europe. Microwave Con$, Madrid, Spain, Sept. 1993, pp. 817-818. “MATLAB: High-performance numeric computation and visualization software,” The Math Works, Inc., Natick, MA. “Microwave SPICE,” HP-EEsof, Inc., Westlake Village, CA. “Touchstone,” HP-EEsof, Inc., Westlake Village, CA. 0. Boric-Lubecke, D. S. Pan, and T. Itoh, “Effect of the increased number of diodes on the performance of oscillators with series-connected tunnel diodes,” in Proc. 6th Inter. Symp. on Space THz Tech., Pasadena, CA, Mar. 1995, pp. 207-215. S. Nogi, J. Lin, and T. Itoh, “Mode analysis and stabilization of a spatial power-combining array with strongly coupled oscillators,” IEEE Trans. Microwave Theory and Tech., vol. 41, pp. 1827-1837, Oct. 1993.

Olga Boric-Lubecke (S’88-M’96) received the B.Sc. degree in electrical engineering from the University of Belgrade in 1989, the M.S. degree in electrical engineering from the Califomia Institute of Technology in 1990, and the Ph.D. degree in electrical engineering from the University of California, Los Angeles, in 1995. She is currently conducting research on miniature submillimeter wave radiometers at the NASA Jet Propulsion Laboratory, under a National Research Council Award. In 1990, she was a research engineer at the Institute for Microwave Technique and Electronics in Belgrade, and in 1991 she was a visiting researcher at the Helsinki University of Technology. Her interests include the characterization of solid-state devices and the development of solid-state high frequency sources.

Dee-Son Pan (M’89) received the B.S. degree in physics from Tsing-Hwa lJniversity, Taiwan, China, in 1971, and the Ph.D. degree in physics from the Califomia Institute of Technology in 1978. In 1977, he joined the Electrical Engineering Department at UCLA as an Assistant Professor, where he is currently an Associate Professor. His research interests are in device modeling, semiconductor physics, and theoretical exploration of new devices. His recent works are in the area of enhancing the power of resonant tunneling diodes as a millimeter and submillimeter source, p-type quantum well device structures, and HBT devices.

Tatsuo Itoh (S’69-M’69-SM’74-F‘82), for a photograph and biography, see the February 1996 issue of this TRANSACTIONS.

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

944

Mode Coupling in Superconducting Parallel Plate Resonator in a Cavity with Outer Conductive Enclosure Feng Gao, Member, IEEE, M. v. Klein, Senior Member, IEEE, Jay Kruse, Member, IEEE, and Milton Feng, Fellow, IEEE

Abstruct- We have carefully studied the mode coupling effect from analysis of the measured microwave scattering parameters of superconducting films using a parallel-plate-resonator technique. Due to its high resolution and simplicity, this technique has been widely employed to identify the quality of high-T, superconducting films by measuring the resonance bandwidth, from which the microwave surface resistance is directly derived. To minimize the radiation loss, the resonator is usually housed in a conductive cavity. Using this method, we observe that a number of strong “cavity” modes due to the test enclosure fall around the lowest TM mode of the superconducting resonator and that a strong interaction between these two types of resonant modes occurs when their eigenfrequencies are close, causing a significant distortion or a strong antiresonance for the resonator mode. To describe this effect, a coupled harmonic-oscillator model is proposed. We suggest that the interaction arises from a phase interference or a linear coupling among the individual oscillators. Our model fits very well the observed Fano-type asymmetric or antiresonant features, and thus can be used to extract the intrinsic ?C. of the superconducting resonator. I.

INTRODUCTION

H

IGH temperature superconductors (HTS) have potential commercial applications in microwave devices [ 11-[7] due to their extremely low high-frequency loss and small dispersion compared with normal metals at temperatures even above the boiling point (77 K) of liquid nitrogen, a readily available coolant that is much less expensive than liquid helium, which is required for conventional superconductors. The electronic properties of HTS can be exploited for use in a variety of high-performance microelectronic components such as analog and logic circuits, vortex flow transistors, multichip module interconnects, delay lines, filters, and infrared detectors. The fabrication of ultrasensitive sensors and production of Josephson microwave mixers with wider electromagnetic wave windows becomes possible with the use of HTS because of the large expected energy gap for these materials. To explore the possibilities of these applications, it is essential to investigate the microelectronic mechanisms and to determine accurately some fundamental parameters of high quality Manuscript received October 19, 1995; revised February 15, 1996. This work was supported by the National Science Foundation under Grant NSFDMR-91-20000, through the Science and Technology Center for Superconductivity under Grant DMR-91-20000. The authors are with the Science and Technology Center for Supercouductivity, Department of Electrical and Computer Engineering, University of Illinois, Urbana, IL 61801 USA. Publisher Item Identifier S 0018-9480(96)03788-X.

superconducting films at microwave frequencies. One key parameter frequently measured is the temperature-dependent complex impedance, Z,(T) = R, j X , , from which the London penetration depth XL ( T ) ,complex conductivity o ( T ) = o1 - j a 2 , superconducting order parameters A ( T ) ,and quasiparticle scattering time T ( T )can be deduced [8]. Microwave surface resistance in principle can be measured in a variety of traditional methods. However, because of the limitation that high-T, materials are highly anisotropic and of the requirements that the HTS thin films can be grown only on substrates such as sapphire, MgO, LaA103, SrTi03, or yttria-stabilized zirconia (YSZ), it is impractical to fabricate an entire resonant cavity from the HTS material for R, measurements. Therefore, one method commonly adopted is to compare the quality factor (Q) of a metal cavity with and without a small superconducting sample enclosed. This perturbation procedure, ideally suited for measuring high R, values, involves a careful and difficult determination of a rather small fractional difference between two large values of separately measured Q’s because the HTS sample contributes only a negligibly small fraction of the total cavity loss. Other methods commonly in use are patterned superconducting stripline [9] or coplanar resonators [lo]. The parallel-plate-resonator (PPR) technique developed by Taber [ 1 1 ] affords an ideal alternative to determine the surface impedance with high sensitivity. It is best suited for measuring R, values ranging from 10 pR to 1 mR around 10 GHz. This method requires no lithographic patterning that may somewhat degrades the material quality. The current distribution over the tested films is quite uniform, making the R, extraction simple and straightforward. Meanwhile, the frequency shift with temperature can be easily identified; hence the penetration depth of the superconducting samples can be inferred [8]. There are several disadvantages, however, with this technique. First, two identical films, both in dimensions and physical properties, are required in the measurement. Second, any metallic residue left over from processing on the back sides or the edges of the substrates can degrade the measured Q . Third, the resonance signal is sometimes difficult to detect because of difficulties in establishing the proper alignment to reach an optimum coupling between the coaxial probes and the resonator. In particular, the parallel-plate mode usually deviates from an ideal Lorentzian shape due to coupling between the resonator and the cavity enclosure, making it

0018-9480/96$05 .00 0 1996 IEEE

+

~

GAO ef al. : MODE COUPLING IN SUPERCONDUCTING PARALLEL PLATE RESONATOR

difficult to determine the intrinsic resonance linewidth. This problem is carefully dealt with m the present work. Because the PPR geometry is in analogy with a stripline configuration where the cavity chamber constitutes the ground plate and the superconducting films form the center strips, the observed cavity modes are associated with both the cavity chamber and the superconducting plates. When two modes (e.g., a resonator mode and a cavity mode) propagate in the chamber, the signal may not reimain in phase and dispersion may occur due to different phase velocities for different modes [ 121, resulting in a considerable signal distortion when their resonant frequencies are close. Although this situation is not desired in practical measurements, it is not always possible to eliminate the undesired mode that falls near the other mode of interest. Therefore, it is necessary to establish an understanding of the coupling mechanism. In this paper, a systematic study of the frequency-dependent resonance spectra of the cavity and resonator is presented. We find that the unusual antiresonance features result from a phase interference of the wave functions between different modes and can be modeled with a theory of linearly coupled harmonic oscillators.

11. THEORY

945

TMo, (E,) Mode

H

Fig. 1. Profile of the field lines for the TMol mode of a parallel plate resonator. For visualization convenience, the distance, d , between the two plates has been greatly exaggerated. In reality, d = 12.5 pm, and a = b = 1 cm, suggesting that the edge effect is confined to a negligibly small fraction of the total area of the PPR.

Note both H and K have no components in the z direction. Equation (3) indicates that the tangential H components near the sample edges are vanishingly small, and that the components of the surface current normal to the edges vanish at the boundaries (ITz= 0 at z = 0; a; K y = 0 at y = 0, b) as expected for open-circuited boundary conditions. In Fig. 1, the field lines for the TMol mode are illustrated.

A. Field Distributions in a Parullel Plate Resonator

B. Determinations of Superconductor Parameters

To understand how the PPR technique works, it is helpful to briefly discuss the electromagnetism of the PPR. The resonator presents open-circuited boundary conditions as illustrated in Fig. 1. Transverse electromagnetic modes can be excited between the superconducting plates to form standing waves. Because the separation (a dielectric spacer of typically d = 10 to 50 pm) between the plates is much smaller than the film dimensions (1 cm x 1 cm) and the microwave wavelength (A 1 cm), three assumptions can be made. First, we can neglect the fringing fields. Second, the electromagnetic fields are essentially constant within the spacer along the normal (or z ) direction. Third, the tangential components of the E field should vanish as required by the boundary conditions for perfectly conducting plates. The solution of the E field distribution is [ 131 n7rx mry E = iEo COS -- COS (1) a b where u and b are the surface dimensions (-1 cm) of the resonator, n and m represent mode indices, and i is chosen normal to the sample surfaces SCI that the modes are named E, or TM,, modes. The magnetic field ( H ) and surface current ( K ) distributions can be obtained [12]:

On substitution of (1) into the Helmholtz wave equation, one can obtain the resonant frequency:

N

where w = 2 n f is the angular frequency and bo is the magnetic permeability of free space. Other terms in ( 2 ) are given by m7r nnx . mxy H , = K y = -jE o cos __ sin POWb b ’ nx nEz mny (3) H y -- - K 2 -’Eo sin - cos -3 a b ‘ p0wa

fo

=

2

/(;)2

+ (y)2

(4)

where vp = c / & is the phase velocity in the medium with a relative dielectric constant F between Iwo plates. In our experiment, E = 2.04, a = b = 1 cm, giving a predicted fo z 10.5 GHz for the TMol or TMlo mode. Note that degeneracy occurs for a = b. In practice, a resonance always has a finite half-power width Af due to microwave losses. When the PPR is housed in a small conductive chamber, the radiation loss due to field leakage out of the PPR increases linearly with the spacing d [ l l ] , [14]. In contrast, the dielectric loss is independent of d and the loss due to R, of the conducting plates varies as l/d. Therefore, if d is small enough, the conductor loss becomes dominant and one can neglect the radiation and dielectric losses. In this case, the unloaded &-factor of the PPR is given by [111 Q=-fo

4f

where WE and W , are the average electric and magnetic energies, respectively, stored in the resonator, and P, is the average power lost in the resonator plates. Therefore, the surface resistance R, can be determined directly from the 3-dB bandwidth Af of the resonance peak. Note that the

946

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

determination of R, is independent of the particular mode excited in the superconducting resonator. The London penetration depth (XL) of the superconducting films can be determine from the resonant frequency. It is the change of XL that results in the frequency shift. In light of the nonzero value of XL, we must account for the effect of field penetration into the superconducting plates by replacing the dielectric constant t of the medium by an effective value [15] (6) where t is the thickness of the superconducting films. Therefore, XL is related to the resonant frequency by

. I -

'

d

t

-----

where TOis an arbitrary temperature. Therefore, once XL(TO) is known, X L ( T )at any temperature can be determined by measuring the resonant frequency f o (T) and then inverting

(7). Finally, the complex conductivity c can be obtained from the surface resistance R, and the penetration depth XL according to 2 s

= Rs

Teflon Dielectric Spacer

+~ P O W O X L

-

(8)

111. EXPERIMENTAL A parallel plate resonator, shown in Fig. 2, was constructed using two unpatterned identical superconducting thin films grown on low loss LaA103 (1 cm x 1 cm) wafers. A thin teflon dielectric spacer ( d = 12.5 pm and t = 2.04) was sandwiched between the superconducting samples to form the PPR structure. The resonator was then positioned at the center of a cavity chamber machined from oxygen free high-conductivity (OFHC) copper. The cavity, a metallic shield meant to minimize the radiation loss, was polished and gold plated in order to reduce surface loss and prevent corrosion. Two nylon dielectric posts (one was spring loaded) were used to press the samples together, insuring that the samples are stationary during measurements and that there was no electrical contact between the resonator and the cavity enclosure. The microwave signal was launched into and was measured at the output of the cavity/resonator assembly using two parallel, semirigid, 50-0 coaxial cables that inserted through the top of the cavity. Both cables were connected to a Hewlett Packard 85 10 microwave vector network analyzer for twoport scattering matrix measurements. Two micrometers were used to adjust independently the distances between the coaxial probes and the resonator to find an optimum coupling. The incident microwave power was chosen at 1 mW to avoid the high-power nonlinear effect [16]. Optimum coupling was obtained in the undercoupled regime so that the loaded and

Fig. 2. Schematic of the parallel plate resonator and the cavity enclosure. A 12.5-pm-thick teflon spacer is sandwiched between two (1 cm x 1 cm) superconducting samples to form the PPR. The microwave swept signal is sent through one port and collected from the other port of the coaxial cables using a HP 8510 microwave vector network analyzer.

unloaded Q were essentially indistinguishable. Such condition usually occurs when the probe tips are within tenths of 1 mm from the top resonator edge positioned between the inner and outer conductors of the coaxial cables. Cryogenic measurements were performed with the setup placed in a liquid helium dewar. A DT-470 silicon diode temperature sensor and a heater resistor were anchored on the exterior surface of the copper cavity. Both the sensor and the heater were connected to a LakeShore DRC-91CA temperature controller so that temperature variation could be monitored and controlled automatically. In addition, there was an intake in the top of the cavity through which helium gas could be sent to purge the test chamber during the cooling and warming processes. This is an important step to avoid water condensation on the inner cavity walls, the probe tips, and the sample surfaces. The frequency-dependent S-parameter data were collected over a large temperature range between 6 K and T, (-90 K) of the YBazCu307 films. I v . RESULTSAND DISCUSSION Because the cavity chamber is not appreciably larger than the PPR, one expects that the test chamber modes would fall close to the PPR modes. The resonator modes are governed by the intrinsic responses of the superconducting samples,

~

GAO et al.: MODE COUPLING IN SUPERCONDUCTING PARALLEL PLATE RESONATOR

1l ..o * h

E --

947

v 1

0.8

h

0.6

Ez 0.2

-

0.3

h;

N

9

0.4 0.2

10.80

10.85

10.90

0.0

9

10

1'1

12

13

14

Frequency (GHz) Fig. 3 . Experimental transmission coefficient IS21I of the test chamber (open circles on every 5th data point), showing three resonance peaks that can be fit well (solid line) by the Lorentz model using an extended form of (12) with three oscillators. 9

whereas the cavity modes are in principle extrinsic to sample properties. However, the presence and the location of the PPR will modify the resonant frequencies of the cavity modes because of the stripline-like geometry.

10

11

12

13

14

15

Frequency (GHz) Fig. 4. Measured microwave resonance for a 400-nm-thick YBa2Cu307 thin-film PPR placed in the cavity chamber at T

(6)

40 denotes the Fourier tr_ansform of the constant potentials

40j 0’ =

1, 2, 3) at the slots, q S is the Fourier domain magnetic potential at the metallization, is the Fourier domain magnetic field normal to the interface (< = H , at y = d ) and %(a,, Pn) is the spectral domain Green’s function of the problem given by coth ymnd 2 7t(a,n, Pn) = > Ymn = 3;. (7)

<

?”n

+

<

For an ac operation, = 0 at the metallization and hence and $s are orthogonal in view of the inner product defined by ( 5 ) . To solve (6), the moment method has been used. The function b;& where Et is a top-hat function expanded is 5 [Le., = defined as (T% in (4), Fig. 2(c)] which is nonzero over the slots only. Like the electrostatic case, due to the orthogonality of the basis , inner product is zero and hence the functions (Et) and T ) ~ their unknown potential does not appear in the system of equations.

<

Computing b, from the system of equations, the normal magnetic field ( E ) at the slots is determined. This field can be used to find the magnetic flux through a slot

where S,(j = 1, 2; 3) is the surface area of the jth slot. The inductance between any two ports can be obtained by isolating the third port [2], [ 3 ] .This condition can be simulated by assigning appropriate magnetic potentials to slots. For example, if 401 = 1, uo2 = -1, and e 0 3 = -1, then the current flows between ports 1 and 3 only. Thus, inductance L13 between these ports is given by 1513 = @1/4 where 4 is the value of current at port 1 (or port 3): I1

=I, = {HedI

-dx a4

=2/ stripwidth (port11

=2J‘

-dz stripwidth (port3)

= 4.

82

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

983

AGROUND

4t&

Port 1

AGROUND

’f

0.0

Freq.(GHz)

60.0

Fig. 5. S-parameters of a T-junction shorted at one port obtained by the equivalent network (eq), by an electromagnetic simulator (em), and by the measurement (exp); (a) the structure (dimensions in pm), (b) the equivalent circuit, and (c) the results. Therefore, the equivalent inductance between planes r1r: and Tsri [L,I L,J, Fig. l(b)] is given by L,I Le3 = ,513 - ElLl - Z ~ L Z where L, is the inductance per unit length of the CPW’s at port p (p = 1, 2, 3). Two other equations for inductances between ports 1 and 2 and ports 2 and 3 can be set up. These three equations yield the equivalent inductances of the area bounded to the reference planes. For the evaluation of L I , L2, and L Y ,the Laplace equation for the magnetic potential is solved for the three CPW lines connecting to the junction.

+

+

111. RE:SULTS The technique presented was initially checked for the convergence. There are three parameters M , AT (the number of Fourier terms), and I (the number of basis functions) which control the values of the components computed. The convergence can be examined by increasing I while keeping 1M and A T sufficiently large. As an example, the convergence of the equivalent inductances of a T-junction with identical arms is presented. In this junction, the

substrate is GaAs ( E ~= 12.9) with 400 p m thickness, metallized at the backside. Slots and strips are 50 and 76 pm wide, respectively, ZI = Zz = 13 = 350pm, a = 2876 p m and J!J = 876 pm. A large “a” is chosen in order to avoid interference of the magnetic wall with the magnetic field of the slot. Slot 3 is discretized into 20 equal subsections in length, slots 1 and 2 each into 10 equal subsections along the length in the z-direction up to the reference planes and into another 10 equal subsections along the length in the 2-direction up to the center strip. All slots are divided into Ii‘ equal subsections in widths. Therefore, the number of basis functions is I = 60Ii. While maintaining M = N = 900, by increasing IC, the number of basis functions was increased and the inductances were computed. The results are shown in Fig. 3. The final values agree with L,1 = 52 pH, Lez = 51 pH, and Le3 = 6.9 pH in Table 111of 141, given for a similar T-junction. Some discrepancies between the two sets of results would be due to a large underneath airbridge in the structure in [4]. To check the accuracy of the technique, the element values of the equivalent circuit of an equal-arm T-junction, Fig. 4(a), were

984

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 44, NO. 6, JUNE 1996

computed. The results are L,1 = Le2 = 19.59 pH, Le3 = 3.57 pH and C, = -5.83 fF. The equivalent circuit was then embedded in a transmission line network, Fig. 4(b), including the air-bridge capacitances (C = 3.5 fF,calculated using the parallel-plate capacitor approximation) and the lengths of the CPW’s from the air-bridges up to the ports defined in Fig. 4(a). Using the MDS software [12], scattering parameters s11, sal, and 3 3 1 of the network were obtained and compared against those computed by a full-wave electromagnetic simulator [5], Fig. 4(c). For the magnitudes of the scattering parameters, there is excellent agreement between the two methods. The agreement for the phases is good at low frequencies, but deteriorates as frequency increases. The maximum phase error is about 25’ at 59 GHz and that has occurred for s11 . The validity of the component values, L,1 = Le2 = 54.90 pH, Les = 1.86 pH, and C, = 9 fF,of the equivalent circuit of a symmetrical T-junction, Fig. 5(a), computed by the present technique was supported experimentally and also confirmed theoretically using the electromagnetic simulator [5], Fig. 5(c). For the measurement, the T-junction had to be shorted at port 1. In this study, the MDS software was used to set up and analyze the equivalent transmission line network of the structure, Fig. 5(b). Some anomalies noted in the experimental results are due to energy leakage into the substrate and resonances and radiation due to long lengths of lines in the structure. In the experiment, the junction was wire-bonded at the reference planes in order to reduce the slot mode effect and as a result, uncertainties arose with the values of the capacitances of these wirebonds with respect to the center strips. Therefore, these capacitances were not taken into account. However, this should not give rise to a significant error in the results, since the wire-bonds are not close to the metallization and hence should have negligible capacitances. Also, a perfect short circuit was considered for the shorted arm.

IV. CONCLUSION A computer technique based on the quasi-static approximation for the determination of the component values of the equivalent circuits of a broad-class of CPW discontinuities was introduced. The technique enjoys the spectral domain formulation in conjunction with the method of moments for approximating the charge distributions at the conductors and the normal component of the magnetic field distribution at the slots. These distributions are used to compute the equivalent capacitances and inductances. The concepts behind the method were illustrated using an example, the CPW T-junction. The convergence of the technique for a T-junction was reported and the equivalent circuit of two other T-junctions were presented. Both the measurements and full-wave electromagnetic simulations supported the accuracy of the results.

REFERENCES [1] M. Y . Frankel, S. Gupta, J. A. Valdmanis, and G. A. Mourou, “Teraherts attenuation and dispersion characteristics of coplanar lines,” IEEE Trans. Microwave Theory Tech., vol. 39, pp. 910-915, June 1991. [2] M. Naghed and I. Wolff, “Equivalent capacitances of coplanar waveguide discontinuities and interdigitated capacitors using a threedimensional finite difference method,” IEEE Trans. Microwave Theory Tech., vol. 38, pp. 1808-1815, Dec. 1990. [3] M. Naghed, M. Rittweger, and I. Wolff, “A new method for the

calculation of the equivalent inductances of coplanar waveguide discontinuities,” in IEEE M T - S Dig., 1991, pp. 747-750. [4] M. Abdo-Tuko, M. Naghed, and I. Wolff, “Novel 18/36 GHz MMIC GaAs FET frequency doublers in CPW-techniquesunder the consider-

ation of the effects of coplanar discontinuity,” IEEE Trans. Microwave Theory Tech., vol. 41, pp. 1307-1315, Aug. 1993. [SI Sonnet Software Inc., 101 Old Cove Rd., Suite 100, Liverpool, NY 13090. [6] P. Silvester and P. Benedek, “Equivalent capacitances of microstrip open circuits,” IEEE Trans. Microwave Theory Tech., vol. MTT-20, pp. 511-516, Aug. 1972. [7] A. F. Thomson and A. Gopinath, “Calculation of microstrip discontinuity inductances,” IEEE Trans. Microwave Theory Tech., vol. MTT-23, pp. 648460, Aug. 1975. [8] M. Riaziat, R. Majidi-Ahy, and I. Feng, “Propagation modes and dispersion characteristics of coplanar-waveguides,” IEEE Trans. Microwave Theory Tech., vol. 38, pp. 245-250, Mar. 1990. [9] W. P. Harokopus and P. B. Katehi, “Radiation loss from open coplanar waveguide discontinuities,” in IEEE MMT-S Dig., 1991, pp. 743-746. [ 101 D. Mirshekar-Syahkal, Spectral Domain Method f o r Microwave Integrared Circuits. New York: Wiley, 1991. [ 111 T. Itoh and A. S. Herbert, “A generalized spectral domain analysis for coupled suspended microstrip lines with tuning septums,” IEEE Trans. Microwave Theory Tech., vol. MTT-26, pp. 820-826, Oct. 1978. [12] Hewlett-Packard Ltd., Cain Rd., Bracknell, Berkshire RG12 1HN, U.K.

High-Power HTS Planar Filters with Novel Back-Side Coupling Zhi-Yuan Shen, Charles Wilker, Philip Pang, and Charles Carter, I11

Absfraci-Novel back-side coupling was used to produce high-power high temperature superconducting (HTS) filters. Several 2.88 GHz, 0.7% equal-ripple bandwidth, 2-pole TEol mode filters were fabricated using T12BazCaCu208 HTS thin films on 20-mil Lakilo3 substrates. The calibrated, measured performance of the filter at 77 K was