Experimental Investigation of Bifurcation Behavior in Nonlinear Microwave Circuits

We present an experimental technique to study bifurcations in periodically forced nonlinear microwave circuits, including even physically unstable (periodic) steady states. The designer specifies a key node in the circuit being studied (often associated with an active device) and the method synthesizes a voltage waveform to match the waveform at the selected node so that no current flows across the interface. This null condition is maintained while a parameter, such as bias voltage, is varied over a specified range. The addition of the external nulling source is able to stabilize a steady state that would be unstable in the original circuit. Various applications are presented.


I. INTRODUCTION
N ONLINEAR microwave circuits under periodic excitation can exhibit a variety of non-standard behaviors including [1]- [3]: multiple steady-states (typically with different stability properties), hysteresis and jump phenomena with respect to variation of a parameter (such as bias voltage) [4], [5], sub-harmonic solutions [6], parametric oscillations, quasi-periodic solutions [7], and finally chaotic solutions [8]. Practicing RF designers have a variety of colorful jargon names for such behaviors, however the modern language of dynamical systems theory [8]- [11] allows for a much more precise description of such phenomena.
Numerical simulation methods such as harmonic halance [3], [12]- [14] are able to capture all of the operation modes in the list above with the exception of chaos. Moreover, a harmonic balance simulation combined with continuation [15] is able to explore different qualitative and quantitative behaviors. A continuation parameter is a real number (typically with a physical interpretation), which is varied systematically over a specified range to drive the device or circuit being studied into different states. Typical examples of a continuation parameter are bias voltage, external temperature, drive level of a signal generator, or even the position of a slug of a mechanical slide tuner. The representation of the steady states versus the continuation parameter, in the form of a bifurcation diagram [8]- [11], should provide a rather complete understanding of all the possible behaviors of a circuit, stability properties, and its boundaries of operation. In the most interesting cases (discussed later in the paper) modification of the continuation parameter drives the circuit * Emecon, LLC † Univerity of Cantabria This work has been funded by the Spanish Government under contract TEC2014-60283-C3-1-R, the European Regional Development Fund (ERDF/FEDER) and the Parliament of Cantabria (12.JP02.64069).
or device under test into qualitatively different statese.g., periodic states with different stability properties.
Numerical simulation combined with parameter continuation is very powerful however in certain practical cases, an accurate simulation model is not available rendering the numerical simulation results anywhere from suspect to useless. The following scenarios illustrate this kind of simulation shortfall.

1) Researchers continue to investigate new semiconductor
devices that operate at microwave frequencies [16]. Typically, the development of an accurate numerical model lags behind experimental results. If a circuit based on a new semiconductor device exhibits some type of bifurcation behavior before a numerical model for it is available, then simulation techniques cannot be applied. 2) Interconnect parasitics tend to be well-characterized on an MMIC [17], however less so on board-level prototypes. Inadequate modeling of such parasitic elements can significantly degrade the accuracy of simulation results, particularly at microwave frequencies. An experimental scheme for studying bifurcation phenomena implicitly includes the effect of such parasitics. 3) Designers are sometimes expected to work with components (particularly at the system-interconnect level) for which accurate modeling data are not available. For legacy components, the data might have been simply discarded or lost. In a more extreme case, the provider of a system-level component might consider detailed modeling data proprietary and not want to reveal them. Bifurcation behavior of such designs can only be investigated experimentally.
These considerations underline the limitations of numerical simulation and the need for an experimental alternative. This paper proposes an approach to parameter continuation that does not require any sort of numerical model of the circuit or device under test. The technique can explore a complete bifurcation diagram for a circuit or device with respect to changes in a continuation parameter, including both stable and unstable states in an experimental setup. The method is noninvasive in the sense that the resulting bifurcation diagram is identical to what would be computed with numerical parameter continuation, assuming an exact simulation model for the circuit were available.
The experimental results can be used, for example, to extract parameters for an approximate simulation model to bring it into closer agreement with measurements and to provide a thorough understanding of the effect of particular circuit parameters on the bifurcation behavior. For a high-power circuit exhibiting bifurcation behavior, the experimental continuation technique can be used to approach a high-amplitude state in a controlled fashion, avoiding a sudden jump in output power that might damage a device under test. In some cases [18], hysteresis behavior is desirable and the method can be used to confirm an implementation. In other cases (in particular, power amplifiers [4], [5] or some recent MEMS devices [19]), hysteresis is undesirable and the method can be used to confirm that it has been eliminated. Other circuits such as microwave frequency dividers have inherently different regions of operation, delimited by bifurcation phenomena. These regions [7], [20] depend on one or more circuit parameters and it is important to identify and quantify them [21]- [23]. Considering such a variety of situations, the method could also have an application in future measurement equipment as it would enable a thorough experimental characterization of the circuit response versus any parameter of interest, including multi-valued regions.
The paper is organized as follows. Section II presents a brief summary of simulation methods for parameter continuation. Section III describes the new parameter continuation method based on a null-current condition and explains how it is theoretically able to stabilize an otherwise unstable periodic steady state. Section IV presents the experimental implementation of this method. Section V demonstrates the method with application to a variety of practical circuit examples. Finally, Section VI summarizes the paper and suggests some directions for possible future work.

CONTINUATION
In order to make the paper self-contained, we present here a very abbreviated description of simulation methods for parameter continuation combined with a steady-state solver. More details are available in the various cited references. We assume the reader is already familiar with a steady-state simulation method such as harmonic balance (HB) [3], [13] in which each periodic steady-state waveform of a circuit is captured with n real numbers (e.g., Fourier coefficients). Figure 1 shows an example of parameter continuation applied to a practical varactor-tuned filter of Figure 12. The circuit is driven by a sinusoidal voltage source at a fixed frequency of f REF = 750 MHz and a fixed input amplitude of +5 dBm. The independent variable in the graphic is the bias voltage applied to the varactor and the dependent variable is the magnitude of the fundamental Fourier term of the voltage waveform across the varactor. For each value of the bias voltage, the simulation has computed a steady-state solution with the specified sinusoidal drive. This circuit exhibits a distinct hysteresis with respect to the bias voltage. In practice, an experimental monotonic sweep up or sweep down of bias voltage would result in sudden jumps in output amplitude at bias values corresponding to the two turning points marked in the graphic (T 1 , T 2 ). Between these two turning points, there is a branch of solutions which are all physically unstable [12]- [14], [24]. With conventional techniques, this branch of solutions can be accessed and studied only in simulation as described in this section. Experimental access to this branch -which we believe to be novel -will be described later in the paper. A different experimental approach to studying bifurcation phenomena has been anticipated recently in the mechanical engineering community. See e.g., [25], [26] for example.

A. General idea of parameter continuation
Keep in mind that the curve shown in figure 1 is actually a projection of a smooth curve in a much higher dimensional space.
The instability of a periodic steady-state solution in the branch between the turning points can be demonstrated in simulation with a conventional transient analysis. If such a steady-state solution is used as the starting point for a transient simulation, one would observe the solution evolve in time by moving away from the initial periodic steady-state and converge to a different periodic steady-state that would be stable under small perturbation.
The instability of the selected periodic steady-state can also be demonstrated with pole-zero identification [6], [27], [28] as shown in Figure 2. The figure presents the evolution of the circuit dominant poles through the right-hand section of the resonance curve in Figure 1, comprised between 2 V and the maximum value 5.8 V. The dominant poles are the ones having the largest real part and, therefore, the strongest impact on the stability properties, stability margins, transient behavior etc. For low values of the first-harmonic voltage amplitude, there are two distinct pairs of complex-conjugate poles (σ 1 ± 2πf 1 , σ 2 ± 2πf 2 ) at frequencies (f 1 and f 2 ) that are incommensurate with the input frequency f REF . At the particular amplitude 3.133 V, one of the pairs of complexconjugate poles splits into two independent pairs of poles at f REF , or, equivalently into two real poles, due to the Fig. 2. Stability analysis as poles move through the right-hand side section of the resonance curve of Figure 1. The analysis is performed with pole-zero identification [6], [27], [28]. The real part of the dominant poles is represented versus the first harmonic amplitude of the node voltage. non-univocal relationship between the Floquet multipliers and the poles [27], [29]. At the particular amplitude of 3.542 V, corresponding to the lower infinite-slope point (T 1 ) of the solution curve, one of the real poles becomes positive and the solution becomes unstable. As the voltage amplitude increases further, this real pole continues to grow, then decreases and takes a zero value at 5.614 V, corresponding to the upper infinite-slope point (T 2 ). For higher voltage amplitudes all the poles have negative real parts and the solution is stable. The section of the curve located between the two points T 1 and T 2 is an unstable region, due to the existence of a positive real pole. The rest of the curve is stable.
Numerical continuation [15] can be viewed as a practical application of the implicit function theorem [30]. Suppose that the complete steady state of a circuit is represented with K = n × m real numbers: x ∈ R K . Here, m represents the number of state variables needed to capture the instantaneous time-domain behavior of the circuit and n = 2N H + 1 is the number of real Fourier components needed to adequately represent a steady-state of each waveform. As a special case, with n = 3, harmonic balance reduces to the describingfunction technique [31]. The mapping f : R K → R K selects the subset of the entire space corresponding to valid states, i.e., x which maps to 0 ∈ R K . Typically f imposes constraints such as Kirchoff voltage and current relationships and I/V characteristics of nonlinear devices.
Suppose furthermore that the circuit is controlled by a continuation parameter µ ∈ [µ 0 , µ 1 ]. A new mapping H : R K+1 → R K is defined such that valid (x, µ) pairs are still those elements that map to 0 ∈ R K . The supporting theory [32] shows that for practical purposes the set of points (x, µ) that map to 0 form a smooth, continuous curve (a one-dimensional manifold) in the ambient space of dimension K + 1.
It is tempting to expect a functional relationship between x and µ. From an engineering standpoint, you turn the µ "knob" and watch the state, x, track different values of µ in a continuous fashion. It is also natural to expect the process to be re-traceable: Move the µ knob in the other direction and the state, x, re-traces its previous path. However, both of these intuitions are completely wrong in the presence of a hysteresis. In order to stay on the zero curve (i.e., those pairs (x, µ) that represent valid states) it may be necessary to adjust both x and µ with either components of x or µ possibly increasing or decreasing locally. In other words, it is not longer valid to think of µ as the independent parameter with x functionally dependent on µ. The zero curve still exists and is continuous, but it exhibits a so-called turning point [15] with respect to µ. The notation y = (x, µ) ∈ R K+1 puts x and µ on equal footing and de-emphasizes the idea that x is a function of µ.
The experimental method presented here will be demonstrated with two different continuation methods: arc-length and a parameter-switching technique, based on the the use of an auxiliary generator [7], [14]. The two methods are described briefly in the following sections. When applied to the varactor filter example, the two methods provide the identical solution curve displayed in Figure 1.

B. Arc-length continuation
An elegant mathematical device is to consider x and µ (hence y) as both depending on arc length r along the zero curve [15], [33], [34]. The starting point of the curve is a point y 0 = (x 0 , µ 0 ) = y(0). Subsequent points along the curve are indicated by increasing the arc length along the curve from this starting point. No one component of the y vector is required to increase monotonically (including µ) but arc length along the curve does increase monotonically by construction.
As mentioned above, the supporting theory shows that the set {y(r)|r ≥ 0} is a smooth curve in the space R K+1 parameterized by arc length r. Moreover for a point y on this curve, the Jacobian matrix ∂H ∂y of dimension K × (K + 1) is full-rank as long as the curve is smooth. The last column of the matrix represents derivatives of the circuit equations with respect to the continuation parameter µ. The Jacobian matrix has a simple and important geometric interpretation. Consider the Jacobian matrix J = J(y) evaluated at a point y on the curve. Since this matrix is full-rank, it is possible [15] to find a unit-length vector v (of dimension K + 1) such that Jv = 0 ∈ R K . This vector is tangent to the zero-curve at the point y; locally it indicates the direction of increasing arc length along the curve.
These considerations can be combined to give a basic algorithm for arc-length continuation: 1) Set µ = µ 0 and find a starting point y = y 0 on the zero curve. Typically µ 0 is chosen to ensure a unique solution outside of any hysteresis region 2) Choose an arc-lengh increment δ.
3) Obtain the Jacobian matrix J evaluated at the point y. 4) Find a unit-length vector v in the null-space of J, for example by performing a QR decomposition on J [35]. 5) Compute a predicted point y advance = y + δv. 6) Consider the hypersphere of radius δ centered on y; obviously y advance is on this hypersphere. If δ is not too large, the zero curve exits the hypersphere at a point y new close to y advance . 7) Use an iterative computation, such as Newton's method, to find the point y new on both the zero curve and the hypersphere employing y advance as a starting point for the iteration.
This basic step is then repeated, advancing in arc length along the zero curve, until the target value µ 1 is reached. In practice, as long as the arc length increment δ is not too aggressive, the algorithm described above will track the zero curve in the form of a series of closely-spaced points along the curve advancing in arc length. Since y advance is close to y new the final iterative step above typically converges in a few iterations and exhibits super-linear convergence. More sophisticated implementations [33] use an efficient adaptive scheme to propose larger or smaller values of δ depending on the nature of the zero curve.

C. Parameter switching with an auxiliary generator (AG)
An alternative simulation technique described in [13], [14] is based on the introduction of an auxiliary generator into the circuit. In HB simulation, this enables the implementation of a parameter-switching method in a two-tier procedure. An early example of this technique appears in [7], [12]. In HB simulation, the designer selects a forcing node (often associated with a key active device) and places a sinusoidal voltage source and ideal transfer function at the AG frequency [13] in shunt from this node to ground as indicated in Figure  3(a). In the analysis of periodic solutions, the AG imposes the first harmonic of the Fourier series of the voltage signal at the node where it is connected. To obtain a hysteresis curve at the drive frequency, such as the one shown in Figure 1, the AG frequency must be equal to that of the driving source. To obtain a sub-harmonic solution it must be a sub-multiple of that frequency. The magnitude and phase of the AG, α AG and φ AG , are unknown variables. A numerical optimization scheme is then used to find the α AG and φ AG that satisfy the non-perturbation condition: where Y (with units of admittance) is the ratio of the AG current to the AG voltage. The continuation method to trace the solution curve versus µ is based on parameter switching. In each section of the solution curve, the variable with the fastest increment (among µ,α AG , and φ AG ) is swept [7], [13], [14] and the other two are calculated to fulfill (1). An important advantage of the AG technique is that it can be applied as a layer on top of an existing commercial simulation tool without a need to change the internal programming. The AG method was applied to simulate the varactor filter described above with the commercial tool Agilent ADS TM . The resulting curve agrees with the one obtained by arc-length continuation shown in Figure 1.
In measurement a different method is going to be used which does not require either an ideal voltage source or an ideal bandpass filter, both of which are difficult to realize in practice. Instead, the method synthesizes a complete waveform f REF

Input
Band Pass transfer function over one period (or multiple periods in the case of a subharmonic response), which must fulfill a null-current condition. Its theoretical foundations will be presented in Section III and its experimental implementation will be presented in Section IV.

III. PARAMETER CONTINUATION BASED ON A NULL-CURRENT CONDITION
As stated above, the continuation method proposed here is to match the complete voltage waveform at the forcing node, not just its fundamental Fourier component, hence avoiding the need for a floating filter which is impossible to implement. Figure 3(b) shows a sketch of a generic circuit inside the dotted rectangle. It is driven by a periodic waveform v REF (t) (typically sinusoidal) of frequency f REF . The circuit is also controlled by a continuation parameter µ ∈ [µ 0 , µ 1 ] represented figuratively as a knob in the figure. Consider a node internal to the circuit with an associated voltage waveform v DU T (t). If the circuit is linear, then this waveform will be sinusoidal at the excitation frequency f REF with period only the magnitude and phase with respect to the excitation source can change. However, in the case of a nonlinear network (especially one containing reactive elements), the situation is more complicated allowing for at least the following possibilities: • v DU T (t) can be periodic, but no longer sinusoidal • there may be different waveforms v DU T (t) if the circuit exhibits multiple steady states (possibly depending on the value of µ) can be a quasi-periodic or chaotic waveform at this node; these cases are outside the scope of this paper.
Suppose it is possible to synthesize a voltage waveform v P (t) equal to v DU T (t) over an entire period at the frequency f REF , or multiple periods in the case of a sub-harmonic response. The substitution theorem of circuit theory [36] then says that the synthesized source can be connected in shunt to the DUT node without disturbing the operation of the circuit. No current flows across the junction because the voltage waveforms are identical. Here we are implicitly assuming that v P is generated by an ideal voltage source. The more realistic case of non-zero source impedance is treated later in this section. The experimental equivalent to the to the nonperturbation condition (1) is that no current flows into or out of the circuit at the forcing node, even though the voltage v P is not identically zero. This is termed the null-current condition. However, even if according to the substitution theorem a solution exists such that v P (t) = v DU T (t) (null-current condition), this solution might be unstable and, therefore, never be obtained in measurement. Actually, in our experimental investigation of bifurcation behavior, we need obtain solutions that are originally unstable, without altering them. Therefore, besides fulfilling the null-current condition, the source v P (t) must be able to stabilize an otherwise unstable periodic steady state. In the following, some mathematical justication for the stabilization mechanism is provided. It is important to keep in mind that a general proof that the experimental method works is not really possible: the success of the method depends critically on the designer's choice of forcing node to which to connect the nulling source. As an extreme case, connecting the nulling probe to ground would not accomplish much. Also, the source impedance associated with v P (t) is important but not critical, as will be shown in a subsequent numerical example. The present implementation uses a source impedance of 50 Ω, which has proven to be adequate for the example circuits presented here. Our experimental experience shows that a large range of source impedance values is usable. On the other hand, a critic of the method could claim that since there is no current flowing through the probe at null, the presence of the nulling probe should have no influence whatever on the circuit or device under test. This is incorrect, as shown by the following analysis. Figure 4 inside the dotted line shows a schematic of the varactor tuned filter presented later as a first practical example. The augmented circuit adds a periodic waveform source v p (t) with source impedance R p . The current flowing through this source is i p (t), which should be zero at a null. Circuit 2 is obviously different from circuit 1 hence it might well have different stability properties. However, for the null-current condition, circuit 1 viewed as a sub-circuit of circuit 2 will have the same steady-state solution waveforms it would have when alone.
Assume the circuit without the nulling source is described by a nodal equation. Note that distributed elements have not been considered here for simplicity, but the extension to distributed elements is straightforward and does not affect the subsequent analysis. where: •X is a vector containing the Fourier components of the m state variables of the circuit •F (X) is a vector of nonlinear resistances and conductances •Q(X ) is a vector of nonlinear charges and fluxes •Ḡ is a vector Fourier coefficients of independent sources • [jω] is the block-diagonal frequency-domain differentiator matrix with blocks of the form jkω with harmonic index k, −N H ≤ k ≤ +N H . The augmented circuit adds a new independent source v p (t) with associated output impedance R p . The waveform of this source is described by a vector of Fourier coefficientsV p , which allows for arbitrary phase shifts with respect to the circuit driving source.
HereĪ p is a vector of Fourier components of the current through the nulling source, which should be zero, andX ′ 1 the Fourier coefficients corresponding to the forced node voltage. The augmented vectorX ′ includes theX from (2) and the additional unknownĪ p .
System (3) can be linearized under a small-amplitude perturbation giving the so-called variational equation: in which the complex frequency s is associated with the perturbation [3]. Here, [jω, s] represents a block-diagonal matrix in which the individual blocks take the form jkω + s where k runs over harmonic indices. Because system (4) is homogeneous, for non-zero perturbation, the augmented characteristic matrix must be singular with determinant zero. The poles defining the stability of the steady-state solution of the augmented system correspond to the values of s that make the determinant zero. Through comparison of system (2) and system (3), it is easy to see that their corresponding variational equations will be different, unless R p tends to infinity. In that situation, the source v p (t) will be actually disconnected from the circuit. Therefore, in the general case, the source v p (t) and its associated resistor R p will have an impact on the stability properties. Note that for the ideal case of R p = 0, the vectorX ′ 1 would agree with V p (in Fourier components), soX ′ 1 would no longer be an unknown of the system.
Here the poles of the augmented circuit will be calculated numerically using pole-zero identification [27], [28]. Such a method was applied to the circuit shown in Figure 4 as implemented in the commercial simulator Agilent ADS TM . First, the simulator was asked to find a steady-state solution in the middle of the unstable fold-back region of Figure  1. The particular solution corresponds to the bias voltage V b = 5.184 V and has the first harmonic magnitude 5.986 V. Then the nulling source was attached with a waveform to makē I p = 0. Poles of the associated Jacobian matrix were tracked numerically as R p was allowed to grow. One would expect that for a sufficiently large value of R p there should be a loss of stability for the steady-state solution. Indeed, as shown in Figure 5, for R p around 5.2 kΩ (since the independent source is no longer able to control the circuit's physical solution), the dominant real pole becomes positive and remains so for larger values of R p indicating unstable behavior. Because in this case instability is associated with a real pole passing through zero and does not involve the generation of any new frequencies, either subharmonic or incommensurate, additional evidence of the stabilization effect of the nulling probe can be provided by an alternative numerical experiment. The first graphic of Figure 1 was generated by a numerical steady-state simulation combined with numerical arc-length continuation. A periodic steady-state solution of a nonlinear circuit containing reactances can be viewed as a DC operating point of a larger circuit which is still nonlinear, but contains no reactive elements. Reference [37] gives an example of such a transformation. Essentially, the circuit is expanded into sampled time intervals over one period of the driving term and a derivative with respect to time for the reactive elements is replaced by a finite-difference approximation implemented with controlled sources. Solution R p = 50 1k 5k 6k 7k (Ω) 6. Sign of Jacobian determinants depending on source resistor value. Stability changes above 5 kΩ are in agreement with the other analysis.
The larger circuit will have a Jacobian matrix associated with the system of equations that defines its DC operating point. Existing theory [38] relates the physical stability of a DC operating point to the determinant of this Jacobian matrix. For the example of Figure 1, three solutions were obtained for the same value of µ representing three periodic steady-states -two stable and one unstable. The Jacobian matrices were obtained for each of these solutions and indeed -following the theory -two have a positive determinant and one has a negative determinant. Observe that the sign of the determinant of a large matrix can be obtained numerically by performing a QR decomposition [35], and then taking the sign of the product of the diagonal elements of the R matrix and the sign of the determinant of the Q matrix.
Next, the circuit was simulated with the addition of a nulling probe with source resistance R p . The waveform of the nulling source was programmed to match the three periodic waveforms at the same node in the three steady-state solutions of the original circuit. Therefore, no current flowed through the augmenting source in steady state. Finally, the Jacobian matrices (K 1 , K 2 , K 3 ) and their determinants were evaluated for the three steady-state solutions of the larger augmented network. As shown in table 6, if R p is sufficiently small, all three solutions to the augmented circuit are stable. However, if R p is larger than about 5 kΩ, the augmenting source is no longer able to stabilize the network, in agreement with the previous analysis.
As has been shown, when increasing R p , the solution forced by the source v p (t) becomes unstable at 5.2 kΩ. Actually, for R p tending to infinity, the circuit should behave exactly as the original one, with the response curve shown in Figure 2. That is, for the particular bias voltage V b = 5.184 V , it should exhibit three coexisting solutions with the respective amplitude values V 1 = 2.662 V, V 2 = 4.986 V (the one that has been stabilized) and V 3 = 5.7 V. Therefore, when increasing R p , the system solutions must evolve from a single one at V 2 to the three original coexisting solutions for R p → ∞.
The evolution of the steady-state solutions versus R p , when forcing the waveform v p (t), has been analyzed using the AG numerical continuation (described in Section II.B). The results are shown in Figure 7, where the magnitude of the node voltage (1) at the fundamental frequency has been traced versus R p . For R p below 1.242 kΩ there is only one solution, agreeing with v p (t). Note that the first harmonic of this solution represented in Figure 7 does not vary with R p since for this solution there is no current flowing through R p for any value. This is why the particular R p value is not critical. However, at R p = 1.242 kΩ, two solutions are generated in a turning point. These are the solutions that will evolve to V 1 and V 3 of the original circuit [not forced by v p (t)] for R p tending to infinity. Observe that these solutions arise only from a significantly large R p value. Fig. 7. Evolution of the steady-state solution of the augmented circuit containing the nulling probe versus its source impedance; As this impedance grows to infinity, the circuit solutions tend to the three solutions (V 1 , V 2 , V 3 ) of the un-augmented circuit.

CONTINUATION USING THE NULL-CURRENT CONDITION
The scheme for experimental continuation based on sensing the null-current condition requires the following components for a circuit driven at a frequency f REF : 1) selection of a suitable forcing node 2) a voltage waveform synthesizer able to construct an arbitrary waveform v P (t) with the same period as the excitation and coherent with it (or, a multiple of this period in the case that a sub-harmonic response is expected) 3) a mechanism to measure the current flow into the forcing node, i N U LL = i N U LL (t) 4) a mechanism to adjust µ under computer control (typically a digital to analog converter) 5) a numerical control algorithm to adjust v P (t) so as to drive i N U LL to zero over an entire excitation period. Assuming these components are available, it is relatively simple to implement a kind of parameter continuation that maintains the null condition for different values of the continuation parameter µ. Observe that there is a numerical computation involved, but it makes no reference to any sort of numerical model of the circuit being studied. The following sub-sections describe (A) the experimental apparatus, (B) the experimental continuation algorithm, (C) the null detector, and (D) Newton's method for getting a starting point.

A. Experimental apparatus
This sub-section lists and describes the components of the experimental scheme shown in Figure 8.
Component #1 is a signal generator (HP8656A) operating at a fixed frequency of f REF above.
Component #2 is a custom-built waveform generator to synthesize v P (t). The signal from the reference generator is split into four channels. As shown in the block diagram, each channel goes through a voltage-controlled analog phase shifter and a voltage-controlled analog attenuator (e.g., Mini-Circuits ZX73-2500). Control voltages for the phase shifters and attenuators are generated by 12-bit digital-toanalog converters (DACs) commanded over GPIB. Channels 2,3 and 4 go through, respectively, a frequency doubler, a frequency tripler and a frequency quadrupler. Each channel is amplified, bandpass filtered and summed in a four-channel combiner. The goal is to get independently controllable terms for a Fourier synthesis. Careful circuit design, shielding, elimination of spurs and attention to power supply filtering were needed to reduce cross-talk. In the final design, the waveform synthesizer can be commanded (over GPIB) to turn on only one tone at its maximum amplitude and all other tones measured at the output of the synthesizer will be lower in power by at least -60 dB. This level of spurious response is consistent with, for example, modern digital AWG implementations such as the Tektronix AWG7000. A portion of the reference signal is tapped off and fed to a comb generator using a step-recovery diode [39] Component #3 is a broad-band, linear power amplifier (Mini-Circuits ZHL-5W) following the summing combiner. The specified third-order intercept point is +45 dBm. At null the amplifier is working into an open circuit. The amplifier must be able to generate enough voltage swing across an open circuit to match the voltage waveform at the probe point inside the DUT. We have not observed any instability with the class-A amplifier specified above. Component #4 is a sampling bridge comprising a broadband dual directional coupler (NARDA 3022), a hybrid combiner, a phase trimmer and an amplitude trimmer. A more detailed description of this component will be given later. But briefly, the coupler measures a forward voltage wave V f at the forcing node and a reflected voltage wave V r . As a point of notation, we write V f (and V r ) to mean a complete periodic voltage waveform represented as a vector of complex Fourier coefficients for each frequency Because of the DC blocking capacitor, the DC term of the Fourier expansion is assumed equal to zero. The VNA is able to measure such coefficients directly in (magnitude, phase) form.
The nulling signal is applied to the DUT through a length of coaxial hardline connected to the output of the dual directional coupler and tip of this hardline used as the reference plane for the measurements of the forward and reflected signals. These two waves are equal at the open-circuited end of a transmission line, so the difference V f −V r can be used as an indication of the null condition. The 180-degree hybrid coupler is used to form the difference of the forward and reflected waves. The difference output of the coupler then represents the null signal. The actual procedure is a bit more difficult and is described in detail in Section IV-C.
Component #5 is a Hewlett-Packard HP8753C vector network analyzer that is able to operate in external reference mode so it can be synchronized with the comb generator output from the waveform synthesizer at different harmonics of the fundamental.
Component #6 is a DAC (HP59501B) which generates µ as a voltage. It has 10 bits of resolution and an output impedance of less than 1 kΩ.

Component #7 is a controlling computer which communicates with the other components over GPIB.
Component #8 is the device under test, in this case the varactor filter.
The use of a dual directional coupler as a null detector is advantageous because such devices are available over a multi-octave frequency range. In general, the experimental continuation apparatus has been designed with a thought towards later extension to higher frequencies.

B. Experimental continuation algorithm
The hardware described above is used to implement the basic continuation algorithm described in Section II but limited to manipulation of a single waveform with four Fourier coefficients. The controlling computer manipulates a vector (x, µ) ∈ R 9 containing four Fourier coefficients in (cos, sin) form and the value of the continuation parameter. The computer uploads the value of µ to the DAC and the eight elements of x to the waveform synthesizer. The algorithm must make only incremental changes in x or µ to ensure that it does not drive the circuit into a qualitatively different mode of operation during the convergence process. A purely numerical implementation of continuation faces the same limitation.
It then interrogates the VNA at four harmonics to obtain the magnitude and phase of the null current as a numerical residual. More precisely, four (magnitude, phase) pairs are obtained from the VNA and then converted into a vector of 8 real numbers representing (cos, sin) pairs; this vector is the actual numerical residual. An estimate of the Jacobian matrix at the present stimulus point is obtained by finite-differencing [35]. The algorithm then applies Newton's method to find the value of (x, µ) needed to achieve a current null. The algorithm then takes a step along the zero curve advancing in arc length (which might not mean increasing µ!) and repeats the process. The present implementation is rather slow because 16 residual measurements, two for each of 8 columns, are needed to get a finite-difference approximation to the Jacobian matrix using central-differences. In principle, it should be possible to measure the four Fourier components of the null current in parallel, giving a significant speedup.
An important practical detail concerns the final amplifier used to apply v P (t) to the device under test. This amplifier sees a wildly varying load as the numerical process converges towards a current null, terminating in an open circuit at convergence. The changing load introduces an additional nonlinearity into the convergence process. This effect can be significantly reduced at the expense of more output power by inserting an attenuator between the amplifier and the device under test. As stated above, we have not observed any amplifier instability with the present experimental setup.
The experimental results to be described in Section V were obtained with the hardware scheme shown in Figure  8. The output of the dual directional coupler (null detector) terminates in a length of 0.141 inch coaxial hardline. This section of hardline is used as a probe to deliver the synthesized waveform to a selected node in the device under test while also allowing the dual coupler to measure the current flow across the interface by calibrating out to the end of the length of coaxial hardline. As indicated in Figure 8(a), the end of this length of hardline is taken as the reference plane for the calibration to be described subSection IV-C. For the cavity resonator of Example 1 to be presented later, the probe was inserted directly into the cavity through a small port in the side. The tip of the probe was used to contact the central conductor of the cavity near the diode. For the microstrip frequency divider example, the probe was connected to the junction of lines L 1 and L 3 , and used a small spring-loaded finger to assure good contact. For the mixer-based divider example, the probe was inserted into a coaxial sleeve near the bandpass filter. Some experimentation, both electronic and mechanical, may be necessary to find a suitable location for, and method to, connect the nulling probe.
For circuits implemented in microstrip, an effective probing technique is to penetrate the substrate from the underside of the board and let the center conductor of the probe touch a point on the microstrip. A good ground connection with a closefitting coaxial sleeve is important. Figure 9 shows the nulling probe contacting the gate of a transistor for the power amplifier example presented later. The access hole should be filled with a mechanically identical but "cold" plug for unprobed operation.

C. Null Detectors and Calibration
For the first attempt at a null detector, a dual directional coupler was connected to a Hewlett-Packard HP54121 stroboscopic sampling scope. This instrument is able to measure periodic signals with a bandwidth in excess of 12 GHz and a dynamic range of approximately 45 dB (8-bit digitization). Following Figure 3, a waveform v P bucks the waveform v DU T . In the case that v P (t) = v DU T (t) over an entire period, no current flows across the interface and the forward and reflected waves are equal. Hence, the null signal was defined as E = V f − V r . Following the footnote above, V f and V r were each obtained as a vector of four complex numbers by reading a time-sampled waveform from the HP54121 over GPIB and performing a Fourier transform. In order to be consistent with notation in the rest of the paper, this procedure can be written is a bit more detail as follows. Signals v f (t) and v r (t) were obtained as (periodic) time-sampled waveforms from the HP54121, then converted into V f and V r by taking the first four terms of a Fourier expansion (ignoring the DC term).
In practice, there is some amplitude imbalance and phase shift between V f and V r introduced by the directional coupler and associated cabling so the actual null signal was defined as E = V f − CV r , where C was a diagonal matrix. Components of C were obtained in a calibration process performed before attempting a null measurement. For the calibration, the nulling probe was open-circuited (actually, "unterminated" might be a better description) and the waveform generator commanded to produce a strong signal at each of f REF , were measured and used to form the ratio given arbitrary V f and V r , the error signal was defined as V f − CV r as above.
An enormous improvement was achieved by using a vector network analyzer as the null detector, along with a passive subtraction bridge shown as components (4) and (5) in Figure  8(a). Specifically, the HP8753C VNA was configured in socalled external reference mode to provide the capability of a modern nonlinear VNA [40]. As shown in the block diagram, a portion of v REF drove a comb generator to make harmonic reference channels at f 0 , 2f 0 , 3f 0 , and 4f 0 . A particular harmonic reference signal was selected by a relay and used as the external frequency reference for the VNA. The very first part of the calibration was to connect the output of the comb generator to channels A and B of the VNA through a power splitter. Then, the waveform source was commanded to generate individual tones near maximum amplitude with the analog phase shifter for each channel programmed to zero volts of drive. The value recorded by the VNA for each of the four Fourier frequencies was stored in the registers of the VNA as "magnitude 1, phase 0". For subsequent measurements, these register values were recalled as appropriate for each Fourier frequency. Note that the VNA is not expected to measure an absolute power or amplitude -everything is relative to the values obtained from the comb generator and the numerical methods simply looks for the smallest possible null. The dual directional coupler was connected as shown to implement a passive subtraction bridge. The amplitude and phase of the two channels were balanced so that the passive subtraction coupler gave a broad-band difference of more than 30 dB between an unterminated and a shorted nulling probe.
As shown in the block diagram of Figure 8(a), a sample of the nulling waveform was applied to channel A of the VNA (called X in the diagram) and the output of the subtraction bridge (called E = E(X) in the diagram) was applied to channel B. Similar to the notation for V f and V r , X and E are used here to denote four-dimensional column vectors of complex numbers as measured by the VNA. The signal E alone was a crude null signal -in other words, driving the magnitude of E to zero indicated an approximate null -but it was possible to improve upon this with a calibration process.
In a calibration phase, the waveform synthesizer was commanded to generate strong individual tones at each of four frequencies and the VNA used to measure the corresponding component of X CAL and E CAL . A diagonal matrix was then defined similar to the above discussion with C k,k = E CAL [k]/X CAL [k], k = 1, ..., 4. The final error signal was defined for arbitrary stimulus as E − CX. This scheme achieved more than 85 dB of dynamic range, defined as the ratio of error signals from an open (i.e., unterminated) probe and a shorted probe. Temperature stabilization of the subtraction bridge (not shown in the photograph) and adequate warm-up time were needed to maintain this high level of performance.
The calibration process just described is obviously related to a classic SOL calibration of a VNA employing a well-modeled OPEN standard. For our purpose, however, there is no attempt to make a quantitative measure of the impedance at the tip of the nulling probe. The numerical method simply tries to drive the null signal as small as it can. The calibration described above should be performed with the same physical probe as will actually be used to connect to the DUT.

D. Starting point for the Newton iteration
Another important consideration for the experimental method is to get a starting point for the Newton iteration at µ 0 (Step 1 of the algorithm outlined in Section II). Newton's method is known to have strongly contractive properties so it is only necessary to get an approximate starting point that is within the so-called basin of attraction of Newton's method [41]. Often the correct phase and amplitude of only the strongest Fourier component are sufficient with all the other components at zero amplitude. In the following discussion, it is assumed that setting µ = µ 0 results in a circuit with a unique periodic steady-state. Hence, if the initial Newton is able to converge it has found the correct starting point for the subsequent continuation study.
One idea is to measure the waveform at the injection node with an active probe. This is demonstrated in more detail in a later section. This can give good results if the node is not excessively high impedance, but again careful calibration of the phase shifts introduced by the active probe must be taken into account. An effective scheme is the following: Measure the waveform at the proposed injection node with the active probe and the circuit under test running at µ 0 and display this waveform on an RF oscilloscope (i.e., a highspeed sampling scope such as the HP54121 mentioned earlier). Note that a spectrum analyzer display is useless because it does not show phase. Then, connect the active probe to the nulling probe exactly at the same point on the nulling probe where it will contact the circuit under test. In other words, implement a kind of RF SPDT switch. Move the active probe back and forth from the forcing node of the DUT to the tip of the hardline used as the reference plane for the nulling probe. Adjust the magnitudes and phases of the synthesized Fourier coefficients of v P (t) to match the waveform on the oscilloscope. Experimental experience shows that this scheme works well as long as the loading imposed by the active probe is not too large.
Another suggestion is to work somewhat indirectly. Measure the output of the circuit under test at its intended output, that is often designed to work into 50 Ω. Use a high-speed RF oscilloscope to get a time-domain waveform. Typically the RF scope would need to be triggered by a sample of v REF . Set µ = µ 0 and connect the nulling probe. Manually adjust the magnitudes and phases of the various Fourier components of v P (t) until the waveform on the output node is close to what it was without the nulling probe. Use the resulting Fourier coefficients as a starting point for Newton's method at µ 0 . This process is somewhat heuristic, but seems to work.

E. Summary of the measurement process
1) Select a forcing node for the Device Under Test.
2) Define a reference phase for the VNA on each Fourier frequency by connecting the output of the comb generator to channels A and B of the VNA through a power splitter. 3) Calculate the four diagonal entries of the C matrix as described in the text. Apply a strong individual tone to the open-ended nulling probe at each Fourier frequency and compute the ratio E/X. 4) Set µ = µ 0 and obtain a reasonable estimate of the waveform at the chosen forcing node. Various approaches to this step have been discussed. 5) Using this estimate get Newton's method to converge to a waveform that accurately drives the null current to zero over an entire period at the driving frequency (or sub-multiple in the case of a sub-harmonic response). In other words, drive the numerical residual vector defined in the text as close to zero in all components as is possible. 6) Begin the arc-length continuation process advancing in arc-length while maintaining the numerical residual vector as close to zero as possible. Stop when the continuation parameter µ ≥ µ 1 (its final value).

INVESTIGATION OF VARIOUS BIFURCATION PHENOMENA
The experimental approach will be demonstrated on a series of example circuits. Presented first is a technique to measure a (periodic) waveform at a very high impedance node in a specific circuit. Simply touching the node with a conventional active probe seriously perturbs the circuit. The null-current technique, however, is able to make an accurate measurement. Next, we present a series of circuits exhibiting two fundamentally different types of bifurcation behavior under periodic drive: hysteresis and sub-harmonic response. We emphasize again that the experimental results shown here have been obtained without any reference to a numerical model and that many of the examples contain a physically unstable region of the corresponding bifurcation diagram that would not be observable experimentally prior to the development of our method. The appendix shows photographs of the actual circuit implementations discussed in this section.

A. Measuring waveforms
The experimental scheme described here can be applied as a general measurement technique for circuits that may or may not display any sort of bifurcation behavior. Suppose it is desired to measure the voltage waveform at a sensitive node in an RF circuit. The standard approach to such a measurement is to use a high-impedance active probe, such as the Hewlett-Packard HP54701A or the Tektronix P6046. However, this class of probe has several practical limitations. The reactive loading imposed by an active probe is on the order of 1 pF which may be excessive for some applications. Moreover, many active probes start to clip for a waveform swing of more than 10 V, preventing their application to RF power circuits.
The use of a nulling probe to infer a voltage waveform is illustrated in Figure 10. This figure shows an attempt to measure the voltage waveform at node 1 in the circuit of Figure  4. This node is very high impedance, at least 5 kΩ at the fundamental frequency. For example, touching this node with a fingertip causes the circuit to completely shut down. Figure  10 shows the output of the circuit across the 50 Ω load R L as measured with an HP54121 sampling scope. The timebase of the sampling scope is set for one period at the frequency f REF (750 MHz in this case). There are actually three traces superimposed here -one trace is for the circuit completely unloaded, one for the circuit loaded by a nulled probe and one for the circuit loaded with the HP54701A active probe. Clearly, the loading imposed by the active probe is perturbing the circuit. In the case where the output waveform is not disturbed, we infer that the measured waveform is correct. The trace of Fig. 10. Output waveforms of cavity resonator influenced by probing a highimpedance node (three waveforms superimposed); The nulled probe can be slid in and out of the cavity only changing the output from trace 1 to trace 2. Figure 11 shows Newton's method converging to a null for a synthesized waveform connected to node 1, the waveform to ||9.940883e-02|| ||6.633292e-02|| ||2.647810e-02|| ||4.777259e-03|| ||1.295721e-03|| ||4.870037e-04|| Newton exits successfully -0.038 -0.682 0.062 0.170 0.041 0.205 0.021 -0.023 Fig. 11. Trace of program output; Upon convergence, the method has calculated Fourier coefficients for the first four harmonics of the nulling waveform ((sin, cos) pairs). There is no DC term because the probe includes a DC block. be measured. The two-norm of the residual (the null current) drops by more than two orders of magnitude (> 46 dB) in a total of 5 iterations. The final program output is a vector of Fourier coefficients ((sin, cos) pairs) for the waveform needed to achieve a current null. The program computes 4 Fourier coefficients requiring 8 real numbers: [a 1 , b 1 , ..., a 4 needed to drive the interface current to null. Note that there is no DC term because the probe signal goes through a DC block.
Upon convergence, it is possible to slide the probe in and out of the cavity -contacting or not contacting the high-impedance node -with almost no observable effect on the output, as shown by the first two traces in Figure 10. The null is not perfect, as shown by the slight difference between traces 1 and 2. Presumably because there is some 5th harmonic energy present in the circuit that is not accounted for in the nulling process.
Finally, a calibration process is used to derive the desired waveform measured directly at the tip of the probe with no loading. Let V f , V r denote the forward and reflected voltages at the probe tip. At null, V f = V r . Let U f , U r denote the voltages measured at the output ports of the dual coupler. For example, an instrument such as the HP54121 can measure an absolute voltage value across 50 Ω. These quantities can be related by a matrix of complex numbers: in which A, B are the main coupling terms and α, β are second-order terms representing the imperfect directivity of the coupler. There will be one such matrix for each Fourier term. After measuring (U f , U r ) for each harmonic term, the corresponding matrix is inverted to give the Fourier expansion of the voltage at the probe tip. Similar to the discussion of the calibration process in Section IV, the calibration performed here is not a complete traditional SOL cal. We make two assumptions: 1) The dual directional coupler is a linear system for each of the four Fourier frequencies with no cross-talk (reasonable unless it is grossly over-driven).
2) The measured waveforms v f (t) and v r (t) (represented by V f and V r ) upon reaching null are identical to what would be measured with the same output from the waveform source, but the tip of the nulling probe disconnected from the DUT (i.e., the condition under which the calibration was performed). If these conditions are satisfied, then a single set of 2 × 2 complex matrices described above (one for each Fourier frequency) are adequate to re-construct the waveform at the reference plane.
Our present implementation is able to synthesize and null four Fourier terms. Higher harmonic terms would be missed by the probe leading to a kind of tracking error. Unlike a conventional active probe however, the nulling probe tends to present the least loading to the strongest harmonic component because the method is able to get a solid null for the stronger terms.
B. Hysteresis behavior 1) Example 1: Varactor loaded resonator: Figure 12 shows a photo of a resonant cavity loaded with a varactor diode. At moderate drive levels (less than about +5 dBm), the device operates as a voltage-tuned bandpass filter, over a range of 700 to 800 MHz with a 3 dB bandwidth of 20 MHz. However, at a larger drive level, the circuit exhibits multiple steadystate solutions as described above. Output coupling is with a small capacitive tap-off. Figure 13 shows an experimental continuation study for this circuit, with the varactor bias as the continuation parameter, obtained with the experimental nulling probe with fixed input drive of +5 dBm. As in the simulation results already discussed, the independent variable in the graphic is the bias voltage and the dependent variable is the magnitude of the fundamental component of the voltage waveform across the varactor.
The experimental scheme was able to trace a complete bifurcation diagram including the unstable branch. The experimental results were taken as a the gold standard and an optimization process performed to modify the circuit parameters of the simulation, achieving reasonable agreement. Exact agreement is not to be expected because of various parasitic reactances that are very difficult to characterize. Suppose the experimental technique is advanced in arc length until it moves into the fold-back section (between the two turning points) and then the probe is disconnected. The circuit will jump to one of the stable periodic steady-states. An early paper by Perlman [42] makes reference to a "bi-stable" region of operation associated with essentially the same circuit used for parametric amplification.
2) Example 2: HEMT power amplifier exhibiting hysteresis: The circuit of Figure 14 is a modification of a practical Lband amplifier presented in [43]. The new design uses a different output network and has been re-tuned to a fundamental frequency of 750 MHz. Gate bias V GS is critical and is introduced via a broad-band bias-T. The modified circuit still exhibits a bifurcation similar to what is reported in the citation. Components C 1 and L 1 form a conventional input matching network. However, as discussed in the cited work, under certain bias conditions, L 1 can form a nonlinear resonance with the input capacitance of the HEMT leading to the bifurcation behavior.  Figure 15 shows experimental results for the circuit (at V GS = −3.77 V) clearly exhibiting hysteresis with respect to input drive level. The results agree qualitatively with the cited work, but an exact comparison is not possible because the circuits are different. Fig. 15. Experimental results for HEMT PA. Output amplitude is traced versus input drive level even through the unstable "fold-back" region. The results agree qualitatively with the cited work.

C. Sub-harmonic behavior
The experimental results in this section were obtained with a periodic drive at the frequency 1500 MHz, but generating output at 750 MHz. The nulling waveform needed a component at both frequencies.
1) Mixer-based frequency divider: Figure 16 shows a socalled regenerative frequency divider. This circuit was used many years ago to perform frequency division before highspeed digital dividers became widely available [44], but has been studied more recently [45] because of desirable phasenoise properties. The input port of the mixer is driven by a sinusoid at the frequency 2f REF . Assume (somehow) that the local oscillator drive is available at the frequency f REF .
The mixer will generate a sum and difference, giving 3f REF and f REF . The filter eliminates the 3f REF image term and the f REF signal is amplified and used regeneratively as the local oscillator drive for the mixer. There must be sufficient input signal at 2f REF to make the loop gain greater than unity, at which point the circuit starts dividing. Figure 16(b) shows an experimental bifurcation diagram for the divider. The independent axis is the drive level at the frequency 2f REF = 1500 MHz. The dependent axis is the magnitude of the output signal at the divided frequency. Observe that there is a sudden turn on behavior, which is typical of a flip bifurcation In this circuit, there is actually hysteresis with respect to the drive level, indicating specifically a sub-critical flip bifurcation [9]. By comparison, example 3 shown next is also a sub-harmonic frequency divider that displays super-critical flip behavior.
2) Frequency divider based on tuned lines: This example is another sub-harmonic divider implemented with microstrip and a varactor diode. With reference to Figure 17 current generated by the varactor diode is not loaded by the input circuitry. Reactances X1 and X2 are chosen to achieve series resonance at 2f REF and f REF , respectively. See [46] for further details of the operation of this circuit. Note that bias circuitry is not shown. Figure 17(b) shows experimental results. The independent variable is drive level as determined by the bias applied to the voltage-controlled attenuator and the dependent variable is the magnitude of the fundamental component of the output signal (at the divided frequency f REF ). Note that this circuit does not turn off as abruptly as the previous implementation, presumably because the quality factor of the microstrip lines is nowhere near as large as that of the resonant cavity used to implement the bandpass filter of Figure 16. In the language of dynamical systems, this example exhibits a super-critical flip bifurcation [9].

FUTURE WORK
The paper has two main conclusions. 1) A periodic voltage waveform applied to an appropriately selected node in a circuit under test and adjusted so that no current flows into the circuit over an entire period (the null-current condition) can be used to stabilize a periodic steady state that is unstable without the addition of the nulled waveform source. However, stabilization is lost if the source impedance of the waveform is too large. 2) Th process described in point 1 has actually been implemented in a working experimental setup and can be applied without any reference to a numerical simulation model. A sensitive measurement of the null-current condition is possible using a passive subtraction bridge connected to a VNA. The application of the method to several practical examples has been discussed. An immediate extension of the ideas presented here would be to move to a higher fundamental frequency, perhaps even milimeter waves. Theoretically, the scheme presented here could be extended to any frequency range in which a multioctave dual directional coupler could be realized. The directsynthesis scheme used to generate the waveform v P (t) should extend to higher frequencies.
For the examples presented here, a single nulling probe is adequate to explore all the steady states of the circuit depending on the value of the continuation parameter. There is absolutely no reason to assume that this is always the case. It is an open theoretical question to characterize circuits for which this is the case or to explore circuits that require more than one nulling probe (e.g., two probes in quadrature) to explore all steady states.
As mentioned above, the speed limitation of the present implementation is the finite-difference approximation to the Jacobian matrix needed for numerical convergence to a current null. Computational experience shows that this matrix is (block) diagonally strong. A matrix element can be interpreted as a ratio of null current to a stimulus voltage, having units of impedance. Is it possible to obtain at least the diagonal elements of the Jacobian matrix using a network analyzer directly, rather than with finite differencing? Can the offdiagonal terms of the Jacobian matrix be interpreted (and measured) in terms of the modern theory of X-parameters [47]?