Analysis and Simulation of the Effects of Distributed Nonlinearities in Microwave Superconducting Devices

26 IEEE TRANSACTIONS ON APPLIED SUPERCONDUCTIVITY, VOL. 15, NO. 1, MARCH 2005 Analysis and Simulation of the Effects of Distributed Nonlinearities in Microwave Superconducting Devices Carlos Collado, Jordi Mateu, and Juan M. O’Callaghan Abstract—This paper presents a comprehensive study of microwave nonlinearities in superconductors, with an emphasis on intermodulation distortion and third-harmonic generation. It contains the analysis of various resonant and nonresonant test devices and its validation using numerical simulations based on harmonic balance (HB). The HB simulations made on test devices show that the closedform equations for intermodulation and third-harmonic generation are only valid at low power levels. The paper also contains examples of application of HB to illustrate that this technique is useful to simulate superconductive devices other than simple test devices, and that the validity of the simulations is not restricted to low drive power levels. Most of the analyses and simulations of this paper are based on electrical parameters that describe the nonlinearities in the superconducting material. These parameters are compatible with many existing models of microwave nonlinearities in superconductors. We discuss the particulars on how to relate these electrical parameters with one of the existing models that postulates that the nonlinear effects are due to a dependence of the penetration depth on the current density in the superconductor. Index Terms—Harmonic balance (HB), intermodulation, nonlinearities, superconductors. I. INTRODUCTION PLANAR high-temperature superconducting (HTS) filters are being actively developed to take advantage of their low volume, reduced insertion losses, and high selectivity. However, due to the dependence of the surface impedance on the applied field, there are still many applications where the use of these materials is limited. For example, intermodulation distortion (IMD) is a serious limitation to the use of HTS components in communication systems. These limitations might be eased if engineers could reliably predict the nonlinear effects of their designs. In this paper, we describe how to make reliable simulations of the microwave nonlinear effects in superconducting devices. The algorithms proposed allow the simulation of complicated devices driven by signals that are not necessarily sinusoidal Manuscript received February 17, 2003; revised January 30, 2004 and April 14, 2004. This paper was recommended by Associate Editor J. E. Mazierska. This work was supported by the Spanish Ministry of Science and Technology through Project MAT2002-04551-C03-03 and by Generalitat de Catalunya (DURSI) through Grant 2001 SGR 0026. C. Collado and J. M. O’Callaghan are with the Universitat Politècnica de Catalunya (UPC), Barcelona 08034, Spain (e-mail: [email protected]; [email protected]). J. Mateu is with the Centre Tecnològic de Telecomunicacions de Catalunya (CTTC), Barcelona 08034, Spain (e-mail: [email protected]). Digital Object Identifier 10.1109/TASC.2005.844134 (such as those used in communication systems).These algorithms use electrical parameters to model the nonlinearities in the surface impedance or in parameters that are closely related to it (like the inductance and resistance per unit length of a transmission line). Extracting these electrical parameters accurately from experiments is a crucial step to obtain good results from the simulations, and this paper analyzes a few simple devices on which this can be done. The paper also contains an initial section (Section II) that relates the electrical parameters used by the simulators with one of the models of nonlinearities in superconductors, which postulates that the nonlinear effects are due to a dependence of the penetration depth on the current density in the superconductor. II. NONLINEARITIES IN SUPERCONDUCTORS Nonlinearities in superconductors give rise to a nonlinear dependence of the electric field on current density. This dependence is analyzed in this section to set the basic equations to be used in the analysis and simulations made throughout the paper. The formulation used is compatible with most physical mechanisms claimed to be responsible for nonlinear effects in superconductors, and allows the prediction of the nonlinear performance of superconductive devices from parameters that characterize the material and are device-independent. A. Nonlinear Conductance and Penetration Depth 1) Intrinsic Dependences: A superconducting material is intrinsically nonlinear due to the dependence of the superfluid density on the current density (where ) [1]. The relative change of on can be characterized by a function defined as (1) which is positive and, in normal operating conditions, much less than one. At moderate current levels, can be approximated by the first nonzero term of its Taylor series expansion, so that it becomes proportional to the square of the current density, i.e., [1] (2) where is the pair-breaking critical current and is a coefficient that depends on the type of gap, the temperature, and the direction of the superfluid flow. In these conditions, there 1051-8223/$20.00 © 2005 IEEE COLLADO et al.: ANALYSIS AND SIMULATION OF THE EFFECTS OF DISTRIBUTED NONLINEARITIES 27 is a dependence of the conductivity and penetration depth on , so that the quasi-particle and superfluid current densities respond to the electric field following the equations: (3) (4) which set the nonlinearity between and the (vector) total current density by which intermodulation products and other nonlinear effects are generated in resonators, filters, and other superconducting devices. The dependence of and on current density can be written as [1], [2] (5) (6) where and the function (which is zero for ) sets how the conductance and penetration depth changes from their small signal values as the current density increases. Reference [1] analyzes the intermodulation products generated by a YBCO microstrip resonator with equations equivalent to (5) and (6), and finds that the intermodulation products due to the nonlinearities in penetration depth [(4) and (6)] are dominant over those due to the nonlinearities in conductance [(3) and (5)]. Furthermore, [3]–[5] show that these intermodulation products are much smaller than those found experimentally, and point out that there have to be other sources of nonlinear behavior, like weak link grain boundaries, with a dominant contribution in their generation. Despite this, the equations used in the intrinsic formulation are the basis of the phenomenological dependences described below. 2) Phenomenological Dependence: There are a number of works [3]–[6] claiming that, as in the intrinsic case, the nonlinear properties of superconductors can be described by a current dependent penetration depth. However, unlike [1], they do not claim that this dependence has to come from a variation of the superfluid density . The phenomenological dependence used in these works would be similar to (6) for (7) where (8) The parameter in (8), instead of being related to the pair-breaking critical current as done in (2), is a scaling parameter that sets the strength of the spurious signals relative to the fundamental one(s), and is obtained by fitting the data from intermodulation or other nonlinear measurements. The values of obtained are significantly lower than what would be expected from (2) , which suggests that, as mentioned in [1], the effects of the intrinsic nonlinearities are not dominant. Reference [7] uses a modulus dependence for (8) instead of a quadratic one, i.e., (9) This dependence can also be inferred from [8] and [9], except that in these works, the data fitting is not done on a penetration depth but on parameters that are specific of a test device, like a current-dependent inductance per unit length in a transmission line. The quadratic and modulus dependences above give rise to a 3 : 1 and 2 : 1 scaling law (on a logarithmic scale) of the power of the third-order intermodulation products versus the power delivered to the device [7], [10]. Other scaling laws are found in experimental measurements [11]–[13], and variations of (8) and (9) might apply to these cases [10]. B. Dependence of Electric Field on Surface Current The nonlinearities described so far relate volume current density with electric field in the superconductor. As in the linear case where surface impedance is used, there are situations where it would be useful to relate the electric field at the surface of the superconductor with a surface current density. This may happen in cavities where superconductors are used as endplates or as material covering the cavity walls. As discussed later in Section IV-F, this also happens in some types of electromagnetic analysis of superconducting planar circuits, where there is no interest in knowing the detail of the current distribution profile along the thickness of the superconducting film. 1) Surface Current Density and Boundary Conditions: In linear media, there are simple relations between surface current density, penetration depth, and volume current density, which have to be reassessed for nonlinear superconductors. To fulfill the boundary conditions at the surface of the superconductor, the magnitude of the surface current density has to be equal to that of the tangential magnetic field at the surface . On the other hand, the volume current density at the surface equals the spatial derivative of the tangential magnetic field in the direction perpendicular to the surface, i.e., The ratio between both current densities has dimensions of distance and can be related to the rate of decay of the magnetic field into the superconductor at the surface [14] (10) In a nonlinear superconducting material in which depends on the local volume current density, the decay of the field in the superconductor might not be exponential, and the value of at the surface does not have to coincide with in general. However, for the particular case of intrinsic nonlinearities in a type II superconductor, [14] shows that . For extrinsic nonlinearities, it is reasonable to assume that, whenever the volume current density at the surface is much smaller than 28 IEEE TRANSACTIONS ON APPLIED SUPERCONDUCTIVITY, VOL. 15, NO. 1, MARCH 2005 (the usual operating condition for most superconducting devices), there is a quasi-exponential decay of the fields in the superconductor and [2]. We will, therefore, assume that the surface current density satisfies (11) and that this fulfills the boundary conditions at the surface of the superconductor, regardless of whether we deal with intrinsic or extrinsic nonlinearities. 2) Time Domain Equations: Equation (11) establishes a one-to-one correspondence between values of volume and surface current density. We can then use as the independent variable to describe the dependence of on the current [i.e., write ] and use (11) to rewrite (4) (with )as (12) Note also that at low currents, in (2), (8), and (9) may be replaced by to introduce as an independent variable in these equations, and define a function that will set the variation of with . In these conditions, using (7) and (12) with instead of , the total electric field at the surface might be split in its linear and nonlinear contributions, i.e., (13) where follows the usual linear equation , and is the additional electric field caused by the nonlinear dependence on the current density (14) where (15) Note that . Otherwise would contribute to the linear part of the electric field. According to the two-fluid model [15], the resistive losses in these equations can be taken into account by including a term proportional to (16) Similarly, a term might be added in (14) to account for the possibility of having nonlinear resistive losses. The simplest choice is to make (17) with , although, as in (16), terms with a second time derivative of the surface current would also give rise to a nonlinear resistive loss. We have chosen to include the term in (17) to be consistent with [16]. In subsequent sections, we discuss how to fit (17) to experimental measurements. For engineering purposes, all that is needed is a determination of and . Comparisons of experimental results with or with published values of may be done through (15), and one of the equations defining [(2), (8), or (9) with ]. 3) Nonlinear Surface Impedance: The concept of nonlinear surface impedance has been used in several works (see, for example [16] and [17]), where surface impedance is assumed to depend on the amplitude of the radio frequency (RF) magnetic field. To relate this concept to (17), we have to note that this equation quantifies the nonlinear electric field and, thus, will account for the departure of from its small signal values. This departure can be quantified by assuming that the surface current is sinusoidal , and finding the ratio between its amplitude and that of the fundamental of the electric field. The real and imaginary parts of will then be given by [11] (18) Therefore, and account for the change of the surface resistance and surface reactance from their values at low RF fields. C. Nonlinear Distributed Parameters in Transmission Lines The nonlinearities in the superconductor will affect the electrical properties of superconducting transmission lines by introducing a dependence of the resistance and inductance per unit length ( and ) on the total current through the line . For planar transmission lines, there is not a general closed-form equation relating and to the nonlinear parameters of the superconducting material. Instead, an iterative numerical procedure based on [18] and [19] is outlined in [1] by which and can be found from the volume nonlinearities in the superconductor described in Section II-A. In this procedure (which we will refer to as WSD), the cross section of the transmission line is divided in small elements, so that the current distribution can be assumed to be constant in each element. The values of , , and the current distribution in the cross section of the line for each value of total current , can then be found by solving the system of coupled-line equations that results from analyzing the interaction among these elements. Once and are known, they can be split in a linear (current independent) term, plus a nonlinear perturbation, i.e., (19) with , . Many microstrip and coplanar transmission lines have been analyzed for various nonlinear functions at small power levels . In all cases, the nonlinearity in resistance COLLADO et al.: ANALYSIS AND SIMULATION OF THE EFFECTS OF DISTRIBUTED NONLINEARITIES 29 Fig. 1. Nonlinear equivalent circuit of a superconducting transmission line segment with length dz . and inductance per unit length and followed the same nonlinear law as . That is, for quadratic nonlinearities (8), we obtain and ( and being constants) and for modulus nonlinearities (9), and (where and are constants). D. Other Models of Nonlinearities in Superconductors In the remainder of the paper, we will use and or and to characterize nonlinearities in test devices and to simulate the nonlinear behavior of superconductive circuits. For engineering purposes, their use is sufficient to take into account nonlinearities in superconductors, i.e., there is no need to rely on the phenomenological models of (7) and the parameters , . There are many other models that can also be expressed in terms of ,,or , , and therefore, are valid for the formulas and simulations presented henceforth. Relevant examples of these models can be found in [16], [17], and [20]–[24]. A very complete listing of these types of models is given in the references contained in [24]. III. IMD AND THIRD-HARMONIC GENERATION IN SIMPLE DEVICES The microwave nonlinear behavior of superconductive devices is usually tested using matched transmission lines or resonant devices such as transmission line resonators, disk resonators, or cavities. Below, we outline the analysis of these devices and give closed-form expressions for the third-order intermodulation and harmonic signals generated in them (Sections III-A and III-B). We also describe experimental data that fits these equations (Section III-C). We will use and to denote the fundamental frequencies, for the intermodulation product at , and for the third harmonic at . Also, the subscripts and will be used to refer to frequency domain variables at and ; and and the subscript will be used in equations that apply to all these frequencies. A. Nonlinear Transmission Lines Fig. 1 shows the nonlinear equivalent circuit of a transmission line segment of length . From this figure, we can derive a nonlinear equivalent of the classical transmission line equations (20) (21) with (22) where is the voltage drop across the nonlinear elements , in the equivalent circuit of Fig. 1. The solution of these equations for the cases of quadratic (, ) and modulus nonlinearities ( , ) is discussed below. 1) Matched Transmission Line: We have analyzed the thirdorder intermodulation and harmonic generation in a transmission line which is terminated with its characteristic impedance at one end and, at the other is fed with two signals at and . We have assumed that the nonlinear effects are sufficiently weak to consider that the fundamental currents do not differ significantly from those obtained in the linear case, thus, we have (23) or, using a phasor form referred to the frequency with and . Additionally, we have made the assumption that the intermodulation products and third-order harmonics are only generated by the action of the fundamental signals on the nonlinear elements and of Fig. 1. This means, for example, that we are not considering higher order effects that could contribute to like the third-order mixing of spurious signals at and . These two assumptions could cause error if the amplitude of the fundamental signals was sufficiently large, or the nonlinearities were sufficiently strong. Therefore, the validity of the expressions obtained is restricted to a maximum power level, as will be apparent in the following sections where we will compare them with numerical simulations that are not subject to these restricting assumptions. Under all these assumptions, and as detailed in the Appendix, the third-order intermodulation current in a low-loss low-disper- sion transmission line with quadratic nonlinearities is (24) where is the characteristic impedance at (see (51) in the Appendix). Similarly, the third harmonics is (25) where is the characteristic impedance at . The amplitudes of and in (24) and (25) are approximately linear with distance as long as and , respectively, which is usually satisfied in planar superconducting transmission lines. Thus, the power of these spurious currents dissipated at the termination should be approximately proportional to the square of the length of the transmission line, as is found in the measurements of third harmonic made in [4], where an expression similar to (25) is derived assuming the restrictions described previously. Equation (25), unlike its counterpart in [4], is not restricted to lines much shorter than the wavelength. 30 IEEE TRANSACTIONS ON APPLIED SUPERCONDUCTIVITY, VOL. 15, NO. 1, MARCH 2005 Fig. 2. Ratio between the power of the intermodulation products and that of the third harmonic as a function of the ratio between reactive and resistive nonlinear terms. Dotted line corresponds to modulo nonlinearities and continuous line to quadratic nonlinearities. The Appendix also details the derivation of the equations for the intermodulation and third-harmonic currents generated in the line with modulus nonlinearities when (26) (27) Equations (24)–(27) allow us to relate measurable quantities like the power of a spurious signal dissipated at the termination of a transmission line with the nonlinear parameters , or , . To do this, one has to start by determining the type of nonlinearities (quadratic or modulus) from the slope of the power of the intermodulation products (or third harmonic) versus the power of the source. For quadratic nonlinearities, the equations above imply a 3 : 1 slope between the power of the spurious signals and the power of the sources when plotted on a logarithmic scale. This is provided that both sources deliver the same power to the line . In the same conditions, a transmission line with modulus nonlinearities gives a 2 : 1 slope. Once the type of nonlinearities is determined, one has to discern between resistive and reactive nonlinearities. This can be done by comparing the power dissipated at the termination at and . For quadratic nonlinearities, this ratio can be found from (24) and (25) (28) where . Likewise, for modulus nonlinearities with (29) where . Fig. 2 plots the power ratios of (28) and (29) as a function of and . From this figure, it is clear that one could determine the resistive and reactive nonlinearities by measuring the ratio of the third harmonic and intermodulation product in a matched nonlinear transmission line. We have not found published experimental data to check these findings since only harmonic measurements are usually done in nonlinear measurements of superconductive matched transmission lines [4], [25], [26]. The validity of (24)–(27) has been confirmed with the numerical simulations described in the following sections. Similar equations describing nonlinear dependences other than quadratic and modulus can be derived following the procedure outlined in the Appendix. Derivation of such equations is necessary whenever the dependence of the power of the spurious signals on the source power does not follow the 3 : 1 or 2 : 1 slopes described previously. 2) Transmission Line Resonator: To analyze the third-order intermodulation products generated along a half-wavelength resonant transmission line of length , we assume that the fundamental frequencies and and the third-order intermodulation product (at ) fall into the bandpass of the resonator. The fundamental current will then be and the intermodulation current will have the same spatial distribution, so its phasor will be of the form , where is the complex value to be obtained following the procedure outlined below. For quadratic nonlinearities, the nonlinear voltage at in a segment of the line can be obtained by substituting into (22) and taking Fourier transforms [as done in the Appendix with (53)] (30) On the other hand, the power generated at along the resonator will be dissipated in the resonator and coupling loads. This dissipation can be calculated through , where is the loaded quality factor and is the average energy stored in the resonator (31) Finally, can be found by substituting and (30) into (31) (32) This equation coincides with an equivalent expression in [1] except for normalization constants, and a factor three in the imaginary part of (32) which, as discussed in [27], comes from different definitions of . The derivation is parallel for modulus nonlinearities; the nonlinear voltage will be (33) where denotes the Fourier transform of at . Assuming that the fundamental signals are identical in magnitude, we obtain (34) COLLADO et al.: ANALYSIS AND SIMULATION OF THE EFFECTS OF DISTRIBUTED NONLINEARITIES 31 The restrictions made in this section are similar to those of the preceding one: Higher order effects that could contribute to are not being considered, and the nonlinear perturbation cannot be too high to affect the amplitudes of the currents at the fundamental frequencies (i.e., there are no compression effects). As mentioned in Section III-A1, this will limit the power at which (32) and (34) hold, and their range of validity will be checked with numerical simulations that are not subject to these restrictions. B. Disk Resonators and Cavities IMD in disk or cavity resonators can be analyzed in terms of the interaction between electric field and surface current density [that is, by (17)] and the functions ,.We will assume that for quadratic nonlinearities , (where , are constants), and for modulus nonlinearities , (where , are constants). The analysis is very similar to the one done in line resonators. In a line resonator, we found the voltage distribution at by means of a Fourier transform of (22) [(52) in the Appendix], and now this has to be done from (17) and its frequency domain counterpart (35) where, as in (33), denotes a Fourier transform. For simplicity of notation, we have omitted the dependence of and with position in the equation above. The next step is to recognize that, at , the spatial distribution of the magnetic field (or surface current density ) will be dictated by the resonant mode, as happened in the case of the transmission line where . All that is needed is to find the amplitude of the magnetic field (the equivalent of ) by balancing the power generated on the surface of the superconductor at with that dissipated in the cavity or coupled out of it . Once this amplitude is known, it is straightforward to determine the intermodulation power coupled out of the circuit. The specifics for transverse magnetic ( ) disk resonators and for cavities with superconducting endplates are given below. 1) Disk Resonator: Superconducting disk resonators are known for their high-power handling capability [28] due to the absence of magnetic fields wrapping around the superconducting film. When a disk resonator of radius is fed with signals at and , the magnetic field on the surface of the superconductor film is , where , and are Bessel functions and [29]. To calculate the electric field generated at , we substitute for in (35). For quadratic nonlinearities, we get (36) The fields at will also follow the distribution of the mode and, therefore, . To find , we can impose an equation equivalent to (31): , where and is the dielectric thickness. The result is (37) Following a similar procedure for modulus nonlinearities (see [7] for a detailed description), and assuming , the result is (38) 2) Dielectric Loaded Cavity With Superconducting Endplates: The analysis and simulation of nonlinearities in cavities with superconducting endplates is also of practical interest since it may provide a way to evaluate nonlinear properties of superconducting films without having to pattern them [10], [13], [30]. An empty cylindrical cavity could in principle be used, but the current densities are too small to produce significant intermodulation fields [31]. Much stronger intermodulation can be expected in dielectric loaded cavities like those used for surface resistance measurements [32], [33]. In these cavities, the transverse electric ( )fields and currents concentrate in the vicinity of the dielectric cylinder reducing the overall losses of the cavity and enhancing the nonlinear effects. In a cylindrical dielectric cavity, the magnitude of the magnetic field on the film surface is equal to that of the surface current density , where is given by [33] (39) In the equation above, is the distance to the axis of the cavity, and are the dielectric and shielding radius, respectively, is the direction propagation constant, and are the direction wave numbers (inside and outside the dielectric), and and are dependent on Bessel and Hankel functions [33]. The spatial distribution of can be found following a parallel procedure to that used for the disk resonator [7], [10] (40) where is the normalized energy stored in the resonator, , , and . Similarly, for modulus nonlinearities with (41) where , , and . C. Superconductor Characterization With IMD Measurements of Resonators Despite the possible advantages of using transmission lines to characterize nonlinearities in superconductors (see Section III- A1), resonators are normally used for this purpose since they can hold high currents with moderate power from the signal sources [11]. All the resonators described above can be used to 32 IEEE TRANSACTIONS ON APPLIED SUPERCONDUCTIVITY, VOL. 15, NO. 1, MARCH 2005 Fig. 3. Power dissipated to the load at the fundamental and intermodulation frequencies for the hairpin resonator described in [34]. Measurements (dotted lines) fit values calculated using (32) and (34) (solid lines). extract nonlinear parameters. These may consist of parameters that characterize nonlinear transmission lines ( , or , ), or parameters that characterize only the nonlinearities of the superconductor . In the examples below, we will describe how to extract these parameters from experimental measurements using (32), (34), (37), (38), (40), and (41). These examples refer to intermodulation measurements which are presumably taken with no compression effects. If that was not the case, (32), (34), (37), (38), (40), and (41) would not be valid, and the power of the intermodulation signals would not have the proper dependence on the power of the sources (3 : 1 slope for quadratic nonlinearities and 2 : 1 for modulus nonlinearities). 1) Experimental Data From a Hairpin Resonator: The nonlinear performance of hairpin microstrip resonators (TBCCO on MgO) is widely analyzed and measured in [34] at 77 K and 7.4 GHz. This reference attributes all intermodulation performance to the inductive term which follows a quadratic-law dependence at low powers, whereas, it is proportional to the modulus of the current at high powers. Thus, we use (32) and (34) to fit the two different regimes, resulting H A m at power levels up to 15 dBm and H A m at higher power. Fig. 3 shows the results of the fit. The values of and that are consistent with the measurements have been obtained by fitting and to the measurements with the WSD method. Our results are A m and A m . Reference [34] only fits the experimental data at low power levels, and gives a value of around A m . This difference with our result can be justified by the fact that [34] defines the nonlinear inductance in (22) out of the temporal derivative, which accounts for a scaling factor of between the result in [34] and ours. As discussed in [34], this value of is lower than the expected one if the nonlinear IMD was only due to intrinsic effects following the model of [1]. On the other hand, the turnover from square-law to modulus nonlinearities is predicted by [14] for intrinsic nonlinearities, but at much lower temperatures than 77 K. 2) Experimental Data From a Disk Resonator: To check the equations for disk resonators, we have used the experimental data in [35], which is taken in a disk resonator of YBCO on LAO at 30 K and 2 GHz. In these data, the slope of the intermodulation products versus input power is 2 : 1 which is consistent with a modulus nonlinear dependence. Thus, we have used (38) to fit the data. Since the magnitude of in (38) depends on , many pairs of values can reproduce the experimental behavior. However, as done in [6], we have assumed that all IMD is due to the reactive term. Under this assumption, the resulting nonlinear parameters are and H m A [7]. This is equivalent to A m by (9) and (15). We note that this is on the same order of magnitude as the result in Section II, even though the measurements are at different temperatures and with different HTS materials. We also note that the found is of the same order of magnitude than the estimated value of the pair-breaking critical current for YBCO [7]. 3) Experimental Data From a Dielectric Loaded Cavity Resonator: We have measured the nonlinear parameters of a state-of-the-art YBCO sample produced by a commercial supplier (Theva) using a dielectric loaded cavity resonator [13]. The YBCO film, 700 nm thick, was deposited onto 10 10 mm MgO substrate with a thickness of 0.5 mm, and was used as delivered, with no further processing on our part. The resonator was similar to the one described in [32], with a rutile cylinder 4 mm in diameter and 3 mm in height. We made IMD power measurements using copper at one of the cavity endplates and YBCO film at the other. As in the disk resonator experiment, the measurement suggests a modulus model. Assuming again that resistive nonlinearities are negligible, we use (41) to obtain H m A . In this case, the characteristic current density is A m , of the same order of magnitude as the value obtained in the disk resonator. Nevertheless, in [13], we show measurements of other samples produced by the same supplier, which show significantly different IMD in the same measurement setup, despite having similar linear properties at low power levels IV. HARMONIC BALANCE FOR SIMULATION OF SUPERCONDUCTING DEVICES There are several important limitations to the analysis and data fitting presented in Section III. A first limitation is the upper power limit at which the equations in Section III are valid. This limitation is set by the higher order effects (like the mixing of higher order spurious signals) not taken into account by these equations. A second restriction is the inconvenience of using the equations of Section III to analyze complicated circuits such as filters with several resonators having multiple couplings among them. Finally, a third fundamental limitation is related to the type of signals used to drive the devices: We have been using sinusoidal signals and, owing to the nonlinearity of the devices, it is difficult to extrapolate the results found so far to other signals that might be more representative of practical applications (like communication systems). For example, [36] finds that the third-order intercept (a parameter derived from two-tone IMD measurements) is not useful to predict the performance of nonlinear devices subject to code-division multiple-access signals. COLLADO et al.: ANALYSIS AND SIMULATION OF THE EFFECTS OF DISTRIBUTED NONLINEARITIES 33 Fig. 4. Equivalent circuit of a microwave device from a discretization (for example, breaking down a transmission lines in short segments). The nonlinear one-ports model the distributed nonlinearities of the superconductor. For the case of discretized lines, the one-ports would be current-dependent resistances and inductances. A numerical simulation technique can overcome these three important limitations and provide an independent way to validate the equations in Section III. In this paper, we describe how to use harmonic balance (HB) to overcome the first two types of limitations. The use of HB to overcome the third limitation (prediction of performance under nonsinusoidal signals) is described in [37]. We show below (Section IV-B) the validation of the equations for IMD and third-harmonic generation in traveling-wave and resonant transmission lines of Section III-A, and describe how to find the limits of validity of these equations. Section IV-C shows how to extend this type of simulation to disk resonators, and Section IV-D describes how to validate and find the limits of the equations given in Section III-B for dielectric-loaded cavities. Sections IV-E and IV-F are related with the second limitation mentioned above, and show how HB can be useful in simulating filters, or planar superconducting circuits with complicated layouts. A. Basics of HB in Superconductors HB is a technique which is widely used in nonlinear simulations of circuits with lumped nonlinear devices [38]. Its effectiveness is due to the independent calculation of linear effects in frequency domain and nonlinear effects in time domain. The use of this technique for the simulation of superconductor nonlinearities requires taking into account its distributed nature, and most of the HB modifications presented here are done to this effect. Otherwise, the method is quite general and is not restricted to specific nonlinear functions, circuit topologies, or power levels. The spatial discretization usually required to handle distributed nonlinearities implies handling an equivalent circuit with a large number of nonlinear elements, and this requires special care in the numerical techniques used to analyze it. We use the term “multiport HB”(MHB) to refer to the versions of HB that are optimized for the large matrices that result when analyzing circuits with many nonlinear elements. In MHB, the linear part of the device can be represented by a linear network whose impedance matrix is known. As shown in Fig. 4, this network is fed with signal sources, and is connected to nonlinear elements. These nonlinear elements might be modeling either a nonlinearity between electric field and current density [as in (17)], or a nonlinearity between voltage and current [as in (22)]. Thus, the steady-state of the device is determined by two vectors and , where contains all the frequency components of the current (or current density) through the nonlinear one-ports, and contains the frequency components of the voltage (or electric field) across them. The components of these two vectors are grouped in subvectors (which we denote as and ) that contain the variables (voltages and currents or electric field and current density) at each port at the frequencies of the sources (and ) and at all frequencies where spurious may be generated by the nonlinearities ( , etc.). These subvectors , are related with the time waveforms and by Fourier transforms. The goal of HB algorithms is to find a pair of vectors , that simultaneously satisfy the equations of the linear network and those of the nonlinear one-ports. The linear equations can be written in matrix form as (42) where accounts for the contribution of the sources to when is a null vector. In other words, since superposition can be applied in linear networks, the effect of the current sources can be taken into account by not connecting the nonlinear one-ports in the circuit of Fig. 4, and calculating the voltage (or electric field) that these sources would produce across these ports . In subsequent calculations, the nonlinear one-ports are connected to the linear network and the sources are disconnected, but their effect is retained in the calculations. On the other hand, the nonlinear one-ports impose a timedomain nonlinear dependence between voltage (or electric field) and current (or current density) through it, i.e., (43) which has to be fulfilled at each port . In other words, (43) results from (17) or (22) after discretizing its spatial dependence. To find the correct values of and , the algorithm follows the following steps. 1) Calculate and store and the impedance matrix at the frequency components of interest (those of the sources, and of the possible spurious signals). 2) Make an initial estimate of by assuming that is a null vector (i.e., negligible voltages or electric fields across the nonlinearities) and inverting (42). 3) Break up in subvectors , and transform them to time domain . The values of across the ports are found using (43). 4) The time waveforms are then transformed back to the frequency domain and a vector is formed regrouping the subvectors . 5) Using and in (42), a new vector is calculated. 6) If is not sufficiently close to the present estimate of ,arefined estimate can be formed by jointly taking into account the present estimate of , and its corresponding (see details below). 34 IEEE TRANSACTIONS ON APPLIED SUPERCONDUCTIVITY, VOL. 15, NO. 1, MARCH 2005 7) Once the new estimate for is calculated, we go back to Point 3. The process ends when no significant differences are found between and . There are several methods of making the estimate of mentioned above [38]. One simple and fast approach [39] is to calculate the estimate for the next iteration as a linear combination of the estimate of and of the current iteration (denoted as and , respectively) using the equation and setting the parameter between zero and one. Slow and stable convergence is achieved for high values of , whereas, fast and potentially unstable convergence rates result for close to zero. A more sophisticated alternative to this is Newton’s method, which has also been tested in [2]. In our simulations, we have found that the simpler approach [39] has a better compromise between computational efficiency and convergence properties. Another important and practical issue is to perform efficiently the changes between frequency domain and time domain. A review of methods to do this can be found in [40]. We have used bidimensional fast Fourier transforms [41] in order to avoid aliasing problems if the fundamental frequencies are very close, as is usual in superconducting resonators. In this way, the computational time is kept independent of the difference between the two fundamental frequencies [42]. B. Nonlinear Transmission Lines In this section, we use MHB to show the validity of (24)–(27), (32), and (34). For this purpose, we model the superconducting transmission lines as a cascade of short segments of length (short in terms of the wavelength), like the one in Fig. 1. In the resulting circuit, we group all its linear elements (the linear resistances, the linear inductances, capacitances, conductances, source impedance, load, coupling network, etc.) to form the linear network of Fig. 4. The nonlinear elements, and , are connected to the corresponding one-ports of the linear network. The results of the MHB simulations for travelling-wave and resonant lines are given below. 1) Matched Transmission Line: We have simulated a microstrip matched transmission line of 0.150-mm width and 47.1-mm length with YBCO ( at 5 GHz) on 0.508-mm-thick lanthanum aluminate ( , ) [1]. The line is fed with two sinusoidal signals of peak source current mA, at GHz and , with MHz. The line is divided in 200 segments per wavelength resulting in 600 segments in total. The nonlinear distributed parameters , are obtained by the WSD method (Section II-C) assuming quadratic nonlinearities with A m [34]. Under these conditions, we get nonlinearities dominated by the reactance [1] and we observe very good agreement between (24) and (25) and MHB simulations. We have also analyzed the case where the contribution of resistive nonlinearities is significant. Fig. 5 shows a particular case where we have considered an artificially high value of resistive nonlinearities (with the value of determined as above) to show that, even in this case, MHB produces good results. Fig. 5 shows the simulated and for quadratic nonlinearities and the corresponding theoretical values from Fig. 5. Current distribution along a matched superconducting line with quadratic nonlinearities for the fundamental, third-order intermodulation product, and third harmonic. Comparisons of the results of (24) and (25) with the numerical simulations. At 5 GHz, the line has a characteristic impedance Z = 52  and an attenuation constant of  =0 : 004 Np 1 m . We have used 1 L =1 : 0775 2 10 H 1 A 1 m . The resistive nonlinearities are taken to be equal to the reactive ones (1 R = ! 1 L ) . (24) and (25). The agreement is very good (discrepancy is less than 1% error at the load) except for a small ripple due to the segmentation. A very close agreement is also obtained between simulations and (26) and (27) for modulus nonlinearities. In this case, the comparison has been made with identical parameters to those of the quadratic case above, except that the WSD procedure has been applied with modulus nonlinearities and A m [7]. As in the quadratic case, we find that inductive nonlinearities are dominant and we make with the value of determined from the WSD procedure. The results obtained are very similar to those in Fig. 1. To delimit the range of validity of the equations in Section III, we have carried out simulations with MHB increasing progressively the input power. Fig. 6 shows the power delivered to the load at the fundamental and intermodulation frequencies assuming quadratic nonlinearities. The inset depicts the error between simulations and (24). Note that above 45 dBm of input power, the change in intermodulation power becomes higher than 10%, and this coincides with the onset of the compression effects in the fundamental signals. This is because, unlike the equations in Section III, MHB takes into account higher order effects that cause compression in the fundamental and spurious signals. 2) Half-Wavelength Resonant Transmission Line: We have simulated a half-wavelength resonator using the same parameters as those of the simulation above. The ends of the line are coupled to a source and a load by means of two 17-fF capacitors, which model two 0.15-mm gaps [43]. The length of the line mm is set to resonate at GHz. The peak source currents are mA (available power dBm), and their frequencies are and with KHz. The source and load impedance are equal to the characteristic impedance of the line . The line is split into segments, although a coarser discretization also gives good results. We have calculated the amplitude of