Analysis and Synthesis of Hysteresis Loops in an Oscillator Frequency Characteristic

A methodology for the analysis and synthesis of multiple hysteresis loops in the frequency characteristic of a voltage-controlled oscillator (VCO) is presented. This is achieved through the coupling of an oscillator inductance to multiple external (passive) resonators with resonant frequencies in the tuning range of the VCO. A possible application to the implementation of a compact chipless radio frequency identification (RFID) system is explored, using the oscillator as a reader and placing the external resonators in the tag. The system takes advantage of the high sensitivity to the tag resonances in the presence of hysteresis, which leads to vertical jumps in frequency versus the tuning voltage. A desired bit pattern would be encoded in the tag by enabling or disabling passive resonances at a sequence of frequencies. In the practical realization, the inductors in the oscillator and the external board are implemented through spiral inductors so that the resonators in the VCO and the tag have strong broadside coupling. The coupling effect is modeled through electromagnetic simulations, from which a linear admittance, representing the coupled subnetwork, is extracted. The multihysteresis oscillator characteristic can also be obtained experimentally through a new methodology able to stabilize the physically unstable sections without altering their steady-state values. Different demodulation methods for reading the tag are discussed.


I. INTRODUCTION
Hysteresis [1]- [4] is often observed in nonlinear circuits (both autonomous and driven) such as voltage-controlled oscillators, injection-locked oscillators and power amplifiers. It is due to the coexistence of stable solutions in a certain parameter interval. In a typical hysteresis (between two stable solutions of the same kind), the solution curve exhibits three coexisting sections and two turning points [1]- [4], as shown in the sketch of Fig. 1(a), where a representative variable X (e.g. a voltage amplitude or the oscillation frequency, in an autonomous circuit) is represented versus the parameter . line) is unstable and, thus, cannot be physically observed. When increasing (decreasing) the parameter from a low (high) value the circuit initially operates in section 1 (section 3) and jumps to section 3 (section 1) at the turning point T 1 (T 2 ). In most cases, the hysteresis is due to an expansion of the circuit negative resistance in a certain interval of the excitation amplitude [5], which enables the fulfilment of the steady-state conditions for three distinct amplitude values. In most previous works, hysteresis was undesired, so several simulation tools [1]- [3] were proposed for its accurate prediction at the design stage.
In harmonic balance (HB) [1]- [3], [6], [7] the circuit variables are represented as a Fourier series, whose coefficients constitute the unknowns of a nonlinear algebraic system, solved numerically. Unlike standard time-domain integration [8], the HB method is insensitive to the physical stability properties of the solution, so with the aid of complementary techniques [1]- [3], it is able to provide the complete solution curve in Fig. 1(a), including the unstable section 2. In the absence of these techniques, when just sweeping the parameter, one obtains either discontinuous jumps or a loss of convergence in the neighborhood of the turning points, where the Jacobian matrix of the HB system becomes singular [1]- [3]. To circumvent this problem, the works [3], [9]- [11] make use of a continuation method, based on the introduction of an auxiliary generator (AG) into the circuit, which enables a simple implementation of a parameter-switching continuation technique [2]. In [12] a mathematical condition, fulfilled at the turning points, is derived, which enables a direct calculation of these points. The works [13]- [14] take a different approach using a contourintersection method, applied to an outer-tier admittance function extracted from HB. On the other hand, the work [15] addresses the experimental characterization of the full hysteresis loop in driven circuits. This requires the stabilization of section 2 [in the sketch of Fig. 1(a)] without altering its steady-state values. The recent work [16] describes a hysteresis phenomenon obtained through the coupling of a free-running oscillator to an external resonator. This behavior had been observed in the so called "dip oscillator" [17]- [18], which exhibits an amplitude change when weakly coupled to a passive resonant circuit, tuned to its free-running frequency. When increasing the coupling, hysteresis is obtained versus the tuning parameter, which, as shown in this work, is due to a multiresonance behavior about the original free-running frequency. As also shown here, this kind of hysteresis is easy to control and several distinct hysteresis loops can be synthesized in the Analysis and Synthesis of Hysteresis Loops in an Oscillator Frequency Characteristic oscillator frequency characteristic. This can be relevant to chipless radiofrequency identification (RFID) [19]- [23] because the vertical jumps [ Fig. 1(a)] give rise to a distinct oscillator response to the external resonators. As another example, in a sensor/actuator, the core oscillator would contain the sensing capacitance and the hysteresis cycles, induced by the coupling to the external resonators, would give rise to an upward or downward frequency jump when this capacitance exceeds or goes below certain levels. The upward/downward jump in a given cycle would provide a control signal to turn on or off a particular switch. This work extends [16] with an in-depth analysis of the oscillator behavior under the effect of several coupled resonators. Multiple hysteresis loops are synthesized in a controlled fashion, for the first time to our knowledge, and the possible application to chipless RFID systems is studied. This implementation will be based on the principle described in [22], where the impedance of the tag influences the frequency of a sweep-tuned oscillator, acting as a reader. Then, the tag bit pattern is detected from the subsequent change to the DC bias current. Unlike the present approach, no hysteresis is involved. In [23] a detailed analysis and simulation of the system in [22] is presented, considering chipless tags with orthogonally-polarized antennas [21]. A practical limitation may be the small impact of the tag resonances on the oscillator frequency. Here the possibility to increase the reading sensitivity by inducing hysteresis, and hence physical jumps, in the oscillator characteristic is explored. The aim will be to use hysteresis loops to encode the bit pattern of the tag. As will be shown, this procedure enables a much more distinct signal for the individual bits encoded by the passive tag than the one resulting from ordinary resonances. The hysteresis loops are synthesized according to the desired bit pattern of the tag, which is encoded by enabling or disabling the resonators that give rise to the distinct hysteresis loops. For instance, Fig. 1(b) and Fig. 1(c) shows the expected responses for "1111" and "1010", where  is the oscillation frequency. Each hysteresis loop occurs in a distinct parameter interval, as shown in the sketch of Fig. 1. Otherwise, in the same parameter interval there could be more than two stable solutions. In the RFID application considered in the paper, this would give rise to uncertainty in the tag reading. Different methods of detecting the encoded bit pattern ("reading" the tag) with different sensitivities and circuit complexities are discussed.
The paper is organized as follows. Section II presents an analytical investigation of the multi-hysteresis mechanism. Section III describes the analysis of the oscillator coupled to the resonators in the external board. Section IV presents the experimental characterization of the hysteresis cycles. Section V discusses the implementation of the tag reader.

A. Oscillator Equation
Let us consider an oscillator loaded with a resonant circuit composed by of an inductor L and a capacitor C osc , as shown in Fig. 2. Now, the oscillator will be placed near a number N of external resonators, each one composed of an inductor, capacitor and resistance L n , C n and R n , where n = 1 to N. It is assumed that the inductors L n can be coupled to the oscillator inductor L. A periodic solution is assumed and the system will be analyzed at the fundamental frequency. The oscillator circuit, excluding the coupled inductor L, will be described in terms of the admittance function Y(V,), as shown in Fig. 2, where V and  are, respectively, the excitation amplitude and frequency. Then, the oscillator is governed by the following equations: where Y T is the total admittance, calculated between the terminals of the inductor (L) in the oscillator circuit, as indicted in Fig. 2, I is the current circulating through L, I n is the current through the inductor (L n ) of each of the external resonators and M n (n = 1 to N) are the coupled inductances, related with the coupling factors k n as n n n M k L L  . The currents I n through the external resonators are calculated from: Note that the coupling among the external resonators themselves has been neglected, which, as demonstrated in Section III, constitutes a reasonable approximation. Each of the equations in (3) can be solved for the particular current, I 1 to I N , in terms of I. This provides: The above expression relates the current through each resonator (L n , C n and R n ) to the induced voltage jM n I.
Replacing (4) for I 1 to I N in the expression for V in (2) and this, in turn, in the main oscillator equation (1), one obtains: where ( ) coup Y  is the input admittance seen from the terminals of the oscillator inductor L (Fig. 2), in the presence of the whole set of coupled resonators: Fig. 2. Sketch of the generic oscillator with an inductance L coupled to multiple external resonators. The exemplary case of a Van der Pol oscillator is considered, with the parameters a = -0.03 A/V, b = 0.01 A/V 2 , R = 200 , C osc = 86 pF. The auxiliary generator used for simulation purposes is also represented.
For a compact expression of Y coup it is possible to define the following impedance terms ( ) n Z  , where n = 1 to N: Then, (5) can be rewritten as: Assuming small values of the coupled inductances M 1 to M n , one can perform a Taylor-series expansion of the second term in (8) about 1/(jL). This provides the following approximate expressions: which can be used to obtain an approximate version of equation (8): At this point, it is interesting to note that 2 ( ) / ( ) n Z L   can be expressed as: Then, (10) can be rewritten as: Thus, under weak coupling, the oscillator behaves as if it were loaded with N series subnetworks, connected in parallel between the coupling node and ground. As also gathered from (11), there is a scaling effect, due to the constant coefficients 2 2 / n M L .
If the resonance frequencies  n are too closely spaced, neighboring resonators may significantly affect one another. In this situation, the Taylor series expansion in (9) may not be accurate enough, due to the large magnitude excursions of the total coupled impedance 1 ( ) ...
In order to avoid this problem, the frequency spacing must be larger than the bandwidth of the resonators in (12), depending on their individual Q factors / n n n L R  . In these conditions, the coupled oscillator can be analyzed by splitting (12) into real and imaginary parts: where the superscripts "r" and "i" respectively indicate real and imaginary part. As gathered from (13) the coupling to the external resonators affects both the real and imaginary parts of the total admittance function Y T .

B. Multi-Resonance Effects
To illustrate the multi-resonance effects, a computational test has been performed, considering the coupling of an inductor L to six external resonators. The inductance values are L = L n = 450 nH. The capacitors C 1 to C N are calculated to obtain equally spaced resonances going from 30 MHz to 40 MHz. The coupling factor is the same for all the resonators, k = k n  n, and given by: k = 0.025. Initially, the resonator resistance is R = R n = 1 . In Fig. 3(a), the red solid curve corresponds to ( )   As shown in the next subsection, to avoid the overlapping of the hysteresis cycles, the maximum of ( ) i coup Y  about each resonance n must be smaller than the minimum of ( ) i coup Y  about the next resonance n+1. This non-overlap condition is violated in Fig. 3(a), since the maximum of the fifth resonance is higher than the minimum of the sixth one. When considering the full oscillator circuit [instead of ( ) i coup Y  only] this will give rise to more than three steady-state solutions coexisting for the same parameter values (overlapping of the hysteresis cycles). The overlapping in Fig. 3(a) is due to the high Q resulting under R = R n = 1 . On the other hand, for excessively low Q, each resonator L n , C n and R n will not give rise to a sufficient difference between the maximum and minimum induced by this resonator on the susceptance ( ) i coup Y  . As a compromise, Fig. 3(b) presents the results obtained for R = 2  (hence lower Q), with clearly distinguished extremes (maximum and minimum) about each resonance frequency n  of the external resonators. Note that there is no overlapping of the full set of maxima and minima.

C. Multiple Hysteresis Cycles
To get insight into the impact of the coupling effects on the oscillator characteristic, only one resonator will be initially considered in (13), with the associated impedance term ( ) n Z  . The real and imaginary parts of the total oscillator admittance (Y T ) are given by: As stated, the coupling effects introduce a scaled version of the admittance of the series network R n , L n , C n (in parallel between the analysis node and ground). For better clarity, the simple Van der Pol oscillator in Fig. 2 will be initially considered. The oscillator admittance ( , ) . Suitable values for the parameters a and b are given in the caption of Fig. 2. In these conditions, system (14) simplifies to: Now, the oscillator will be tuned by varying C osc . It is assumed that for the particular capacitor value C osc,n , and in the absence of coupling effects (M n = 0), the oscillation frequency is 1/ Y  can also be zero at two additional frequencies, below and above n  , as easily gathered from (15). Replacing the three frequencies in ( , ) r T Y V  one obtains three different amplitude values, one for each frequency, which explains the hysteresis induced by linear coupling effects. Note that, in the general case (14), the amplitude dependence of the imaginary part ( , ) i T Y V  , will give rise to a small deviation from the original free running frequency  n . Now, the Van der Pol oscillator will be analyzed in the presence of N coupled resonators. This is done by replacing for distinct values of C osc (acting as the tuning parameter) one obtains the curves in Fig. 4  In the presence of the six external resonators, and under variations of C osc , one obtains six hysteresis cycles. This can be seen in Fig. 5. For good accuracy the curve has been obtained in commercial harmonic balance (Keysight Advanced Design System) [1]- [3], [6], [7] with a number (NH) of harmonic terms NH = 7, introducing an auxiliary generator (AG) into the circuit [3], [9]- [11], as shown in Fig. 2. The AG is an independent voltage source in series with an ideal bandpass filter, operating at the oscillation frequency  AG =  and amplitude V AG (both must be determined in the analysis procedure). The AG is connected in parallel at a sensitive circuit node (Fig. 2), such as a device terminal, and must fulfill a non-perturbation condition [2]- [3], given by the zero value of the ratio between the AG current and voltage Y AG = I AG /V AG = 0. In commercial harmonic balance [3], [9], the condition Y AG = 0 is achieved through optimization of the AG frequency and amplitude, through it is also possible to optimize other quantities.
The AG allows passing through turning points, which cannot be done through a simple parameter sweep, as indicated in the introduction. This is achieved by switching the analysis parameter to an AG variable (either frequency or amplitude) in the multivalued interval. In the case of Fig. 5, the AG frequency has been swept (which agrees with the oscillation frequency), solving Y AG = 0 in terms of the oscillation amplitude, V AG , and the physical parameter C osc . With this parameter switch, the infinite-slope points become zero-slope points, without any simulation difficulty.
To emphasize the dependence of the hysteresis cycles on the resonance effects induced by the coupling, the value of the capacitor C osc has been represented versus the AG frequency. This can be compared with the imaginary part of the total admittance, ( ) i T Y  , represented on the right axis. In the calculation of the oscillator curve, three different values of the coupling factor have been considered: k = 0.02, 0.04, 0.06. As gathered from the figure, the six multi-valued regions are located about the resonance frequencies of the six external resonators. The excursions due to the resonance minima and maxima increase with the coupling factor and will give rise to overlapping for too high k. Fig. 5(b) shows the oscillation amplitude versus the tuning capacitor C osc for the same three k values.
To summarize, hysteresis arises because of the particular form of variation of the susceptance about the original oscillation frequency  n under weak coupling conditions (see Fig. 4). If the active device can supply negative conductance at the three resonance frequencies, the steady-state oscillation conditions, Re(Y T ) = 0 and Im(Y T ) = 0, will be fulfilled at three distinct solution points. This will generally be the case, since the device will exhibit negative conductance in a certain frequency band and the two additional resonant frequencies [at which Im(Y T ) = 0] are close to the original one. The condition Re(Y T ) = 0 will be fulfilled for a different amplitude in each case, since, due to the coupling, Re(Y T ) exhibits a frequency dependence.
Regarding the stability properties, as shown in [3], the oscillator exhibits a dominant real pole proportional to det and [JY T ] is the Jacobian matrix of the admittance function Y T with respect to the two state variables V and . Expression (16) formally agrees with the one derived in [24] under a fundamental-frequency analysis. Thus, the solution will exhibit a negative dominant real pole for det[ ] 0 T JY  , corresponding to the stability condition derived in [24]. In general, the second term of det[ ] T JY in (16) is much lower than the first term. In most cases, the standalone oscillator will fulfill: First condition implies a reduction of negative conductance with the excitation amplitude, as occurs in all physical devices from a certain amplitude value, and second condition implies an increase of the susceptance with frequency as in a Foster network. Under weak coupling to an external resonator, one will have three steady-state solutions, as explained above. Due to the form of variation of the susceptance (Fig. 4)   at the original oscillation. Because the coupled resonator is linear and the coupling is weak, this coupling will have only a small effect on the total conductance. Thus, the coupling will destabilize the original oscillation and will give rise to two stable oscillations at a higher and lower frequency,  The edges of the multi-valued regions of the oscillator solution curve are determined by a turning-point condition [3], [12] at which the Jacobian matrix of the oscillator equations becomes singular. In terms of the total admittance function, this singularity condition is given by: As gathered from the previous study, the sign of the derivative / i T Y    will undergo two changes about each resonant frequency  n . It will change from positive to negative and then to positive again (Fig. 4). As stated, in ordinary oscillators, the second term in the determinant of (17) is smaller than the first term, so the change of sign of /

III. SYNTHESIS OF THE OSCILLATOR MAGNETICALLY COUPLED TO THE RESONATORS IN AN EXTERNAL BOARD
In the proof-of-concept implementation, the Colpitts oscillator in Fig. 6, implemented with a bipolar transistor, has been considered. Varactor diodes are used for tuning. The choice of a Colpitts circuit using a bipolar device is convenient for the VHF frequency range (40-50 MHz) of the current experimental setup, but other oscillator topologies could be used as well. The goal will be to obtain four hysteresis cycles in the oscillator tuning range between 40 MHz and 50 MHz. Note, that as reported in [16], the same hysteresis effect has been observed experimentally at a center frequency in excess of 600 MHz. The circuit of Fig. 6 could be extended to sweep over a larger frequency range (e.g., to allow more bits in the tag) with a band-switching arrangement.
As shown in the following, the coupled subnetwork (composed by the oscillator inductor L and the external resonators) can be implemented separately from the oscillator core. This independent design is possible because of the relatively small coupling effects required for a multi-hysteresis curve with no overlapped hysteresis cycles.

A. Implementation of the Resonators
The coupled resonators will be synthesized using planar spiral inductors and lumped capacitors. The planar nature of the inductors is well suited for an RFID tag implementation. Their moderate quality factor will help prevent an undesired overlapping of the maxima and minima of ( ) The spiral geometries required to implement the oscillator inductor L and the inductors L t in the external board have been initially estimated from the inductance expression based on current sheet approximation [25]- [26]. This depends on the spiral number of turns, average diameter, fill ratio and some constant coefficients, determined by the spiral shape (square, hexagonal, circle…). As stated, four resonances will be synthesized in the range 40 MHz to 50 MHz.
The external inductor value is set to L t = 474 nH. Fig. 7(a) shows the symmetrical layout of a board with four squaredshaped spiral inductors. These board inductors have n T =7 turns and external and internal diameters d out = 17.5 mm and d in = 1 mm, respectively.
The oscillator inductance is L= 469 nH. Although this value is similar to the previous one, it has to be implemented in a different manner to enable near-identical coupling effects in all the inductors of the external board. This is done by reducing the number of turns to n T = 2 and increasing the external and internal diameters to d out = 59.25 mm and d in = 52 mm. The whole coupled structure is shown in Fig. 7(b), where the oscillator spiral inductor, with a square shape too, placed over the board inductors, can be noted.
As gathered from the inspection of Fig. 7, the central orthogonal axis of the spirals L t in the external board exhibits a lateral displacement x = 14.7 mm with respect to that of the spiral L. This will give rise to a reduction in the coupling effects that can be compensated, if needed, by approaching the inductor L to the external board. The coupling coefficient k of each inductor can be easily estimated following the expressions in [27]. For x =14.7 mm and a vertical distance of 1.5 cm, the estimated coupling coefficient k is 0.02, in the order of those considered in the resonance analysis of Section II.

B. Electromagnetic Simulation of the Coupling Effects
For an accurate calculation of the coupling admittance, an electromagnetic simulation [28] of the coupling of the four spirals in the board and the big-sized spiral is carried out. Initially, a 5×5 scattering matrix describing the whole passive configuration is obtained, defining Port 1 between the terminals of the oscillator spiral and Port 2 to Port 5 at the locations where the capacitors should be connected (at a later stage) to the board spirals. For a distance d = 22.5 mm between the spiral L and the board, the magnitudes of the scattering parameters S j1 , where j = 2 to 5, are shown in Fig. 8, where they have been represented versus frequency. As can be seen, the power transfer is quite similar for all the spirals.   . Imaginary part of the input admittance Y coup () seen when defining a single port between the terminals corresponding to Port 5 in Fig. 7(a). The results of the electromagnetic simulation are compared with those obtained with the theoretical values of L and L t , and the estimated coupling factor k (in dashed line).

C. Oscillator Solution
The coupled oscillator circuit will be analyzed with an original method, which should facilitate as many electromagnetic tests as needed, without having to repeat demanding harmonic-balance simulations. This is because the total admittance function is calculated as: where, as shown in Fig. 6  is obtained connecting an auxiliary generator [13], [14] in parallel at the node RN, as shown in Fig. 6. The high value parallel resistor R ∞ is used to prevent any convergence problems. A triple sweep is carried out in V tune ,  and V, performing a harmonic-balance simulation at each sweep step, with as many harmonic components as desired. At each steady-state solution, both the real and imaginary parts of Y T must be equal to zero, which, for each V tune , provides two real equations in two unknowns V and . The solutions of this equation system can be obtained through the following geometrical procedure.
For each V tune , the function Y T in (18) is calculated versus the excitation amplitude V and frequency . The real and imaginary parts of Y T provide two surfaces in the spaces V, , Re(Y T ) and V, , Im(Y T ), respectively [14]. As an example, Fig. 10(a) and Fig. 10(b) show the two surfaces obtained in the Colpitts oscillator for V tune = 4.5 V, when considering a distance d = 22.5 mm between the oscillator and the board.
For each V tune , the solutions ( ) tune S V of the complex equation Y T = 0 are the intersections of the two contours in (19) and (20), that is, ( The intersections obtained for V tune = 4.5 V are shown in Fig. 10(c) and Fig. 10(d). The latter presents an expanded view about the region where the two contours intersect.
Note that the method will be valid provided the oscillator waveform at the resonator node (RN) has a limited harmonic content, since it neglects the coupling effects at higher harmonic terms. This condition is valid for reasonably high-Q inductor. The same assumption is made in the experimental implementation of Section IV. intersect at three points (as in the case of V tune = 4.5 V), which provides three steady-state solutions. In Fig. 11(b), the solution points obtained through the contour intersection of Fig. 11(a) have been represented in terms of  versus V tune . In Fig. 11(b), the solution curve is compared with the one obtained through a HB analysis of the whole coupled system (with NH = 7 harmonic terms), with excellent results. The contour-intersection method enables a straightforward calculation of the oscillator solution curve under variations in any parameter affecting the electromagnetic simulation, since there is no need to recalculate ( , , ) tune Y V V  . Comparing Fig. 9 and Fig. 11(b), one obtains a good agreement between the resonances in Y coup () and the three hysteresis cycles in Fig. 11(b). Actually, the observation of the hysteresis cycles is quite independent of the particular frequency characteristic of the oscillator circuit, due to the small coupling effects.

A. Experimental apparatus
The experimental results presented here have been obtained from an updated version of the technique described in [15]. Specifically: 1. The previous work treated only the case in which the device under test (DUT) was a driven system; here, we have extended the idea to the case of an autonomous DUT (i.e., an oscillator).
2. An entirely new numerical method of path-following in n-dimensional space has been implemented using simplicial decomposition [29], [30] which does not require estimating a matrix of partial derivatives, hence is more suited to experimental work. A description of the implementation of a simplicial decomposition algorithm is provided in an appendix. Fig. 12 shows the experimental apparatus and its block diagram, with the Colpitts oscillator (inside the dotted lines) implemented with a bipolar transistor as the autonomous DUT; the tuning voltage V tune is generated by a digital-toanalog converter, DAC1. A signal is injected into the emitter of the Colpitts oscillator from an external VCO whose frequency is controlled by DAC3. The output of the injection source is buffered to prevent it being pulled by the DUT. This is an important practical point. The (buffered) output of the reference VCO is applied to the injection port of the DUT through a voltage-controlled attenuator driven by DAC2. Moreover, the injection current goes through a floating sense resistor (R sense ) connected to a differential active scope probe with a 50-Ohm output. A portion of the reference VCO output is applied to the reference port of an HP8753E network analyzer configured in "external source" mode. Finally, the output of the differential scope probe is applied to the A channel of the VNA, which functions as a phase-sensitive detector with large dynamic range.  Fig. 12. Experimental apparatus and system diagram; Three DACs are controlled by a computer which also interrogates the VNA over GPIB. The injection current (ideally zero) is measured across a series sense resistor with a high-impedance differential scope probe. The VNA serving as a sensitive null detector tracks the frequency of the reference VCO.

B. Numerical method
The entire measurement setup implements a mapping f from R 3 into R 2 . The input of this mapping is (V tune , V ref , V ampl ), which are, respectively, the tuning voltage of the DUT oscillator, the tuning voltage of the reference VCO and the control voltage for the amplitude modulator (a voltage controlled attenuator). The amplitude of the reference oscillator is fixed but the control voltage V ampl adjusts the amplitude of the injection signal to match the voltage waveform of the DUT. The output is Re(I inj ), Im(I inj ), where I inj is the complex current circulating through the resistor R sense . The mapping 3 2 : can be compactly expressed as: According to standard theory [29], the zero-set of this mappingi.e., those triples that drive the real and imaginary parts of the injection current to zeroforms a one-dimensional manifold (a curve) in the threedimensional space. The curve-tracking process starts with initial values for V tune , the tuning voltage of the reference VCO, V ref , and the amplitude setting of the reference VCO output V ampl to achieve a nulled current (in both real and imaginary parts); call this point 0 x . Then, a numerical algorithm (to be described) advances along the zero-curve taking small steps in arc-length along the curve. In other words, arc-length is the quantity that advances monotonically during the tracing process. Particularly in the case of a foldback or hysteresis response, any component of can go up or down during the trace, but arc-length increases along the zero curve by construction. As described in the previous work [15], an application of an injection signal with nulled current does not change the state of the DUT, but can change its stability properties. In the case of an autonomous DUT, the frequency and phase of the oscillation of the DUT are determined by the reference oscillator (VCO) because the DUT injection-locks to it. In this sense, tracing a curve in the autonomous case is easier than in the driven case, because the phase of the injection signal does not need to be modulated during the tracing process. Fortunately, the signal at the emitter of the Colpitts (Fig. 12) is close to sinusoidal because of the high-Q tank inductor allowing the use of a sinusoidal VCO as the reference source. A more advanced implementation would attempt to synthesize a VCO waveform that is an exact match to the emitter waveform, possibly including higher harmonics. An alternative detector for the nulled injection current could probably be realized at VHF with a digitizing oscilloscope. The digital scope would be triggered from the VCO and the digitized waveform passed through a Fourier transform to extract real and imaginary parts at various harmonics.
A high-level description of an algorithm for tracing such zero-curves using simplicial decomposition is provided in Appendix A.

C. Experimental results
In Fig. 13, the simulation results are compared with the experimental measurements obtained when sweeping over the V tune range 0-10 V for a frequency sweep 40-50 MHz. The frequency exhibits hysteresis jumps because of the coupling to the passive tag resonators. The controlling computer sweeps monotonically V tune in small steps up and then back down. At each step, the frequency of the Colpitts oscillator is monitored with a 500  probe on the collector of the transistor driving a counter with a General Purpose Interface Bus (GPIB). Qualitatively similar results can be obtained with a Colpitts circuit in which the tuning is accomplished by varying the value of a mechanical capacitor (which is linear at any fixed capacitor value) demonstrating that the hysteresis is not a nonlinear effect introduced by the varactor. Fig. 14 shows experimental results using passive tags with different bit patterns, "1111" in (a), "0101" in (b) and "1101" in (c). The results are compared with simulations in the three cases. In Fig. 14(a) the same sweep used in Fig. 14(b) was repeated with the entire stabilization mechanism of Fig. 12, using the arc-length continuation algorithm. The resulting traces are shown in red and can, indeed, traverse the unstable portions of the hysteresis jumps.
As shown in Fig. 15, the measured stabilized and unstabilized traces do not agree exactly (i.e., even in the physically stable portions of the curve). We attribute this discrepancy to a parasitic effect at the injection port which has not been completely nulled out. Another possibility concerns the voltage waveform at the emitter during the open-loop sweep. This waveform has clearly some second-harmonic content, approximately -17 dBc. As described above, a more sophisticated implementation would program an injection source to match both the fundamental and second harmonic (hence, nulling them) along the zero-curve. We are presently investigating such an improved implementation. V. READING THE TAG Several methods of reading the tag have been considered and implemented. A simple technique would apply a stepped voltage ramp waveform as the tuning voltage of the VCO and measure the oscillator frequency at each voltage step. The reader would ramp up and down over the same tuning range and look for an obvious frequency difference at the frequencies assigned to bit positions of the tag for the ramp-up versus ramp-down sweep. A significant frequency difference would encode a "1" bit and little or no frequency shift would encode a "0" bit. The ramp-up and ramp-down frequency measurements for a "0" bit might be slightly different due to noise, but nowhere near as large as when the passive tag resonator induces hysteresis at the bit frequency. The experimental results presented here have been obtained with such a scheme using an HP5386 counter connected to a computer over GPIB. This scheme seems implementable in a portable tag reader, but would require a micro-processor control and a reasonably accurate frequency counter operating at the VCO frequency. This might be too expensive for some applications.   We have implemented an alternative technique employing dual oscillators as shown in Fig. 16. The tuning voltages of two VCOs are driven from the same sawtooth waveform. The oscillators track but with a deliberate offset in frequency as shown in Fig. 17. The offset is necessary to prevent the two oscillators to injection-lock to each other when coupled through the mixer.
Oscillator 1 in Fig. 16 will be coupled to the tag resonators while oscillator 2 is designed to have very little coupling. For example oscillator 2 could use a toroidal inductor which is self-shielding. A mixer and low-pass filter are then used to extract the difference in frequency between the two oscillators as a function of the common tuning voltage. In the absence of a tag, this frequency difference is close to a constant value. Our present implementation uses an analog linearizer network to achieve some improvement in linearity between the tuning voltage and frequency of the VCO but it is a very simple circuit. A more sophisticated implementation should be able to get the two traces of Fig. 17 very nearly straight and parallel. However, when a tag is coupled to oscillator 1 there are distinct jumps in the frequency difference corresponding to "1" bits in the tag. Fig. 18 shows a read of the bit pattern "0101" at distances of 1.5 cm and 2.5 cm. The present implementation measures the frequency difference with a counter but --considering that the difference frequency is small compared to the tuning frequencies -it should be possible to use an analog frequency-to-voltage converter [31] to generate a DC voltage proportional to the frequency difference for each tuning voltage. Such circuits are simple and have a fast response allowing for a small period to the sawtooth signal. Table I presents our implementation of the RFID system is based on the principle described in [22], where the impedance of the tag influences the frequency of a sweep-tuned oscillator, acting as a reader.
Finally, a more sophisticated method would actually perform a stabilized sweep with a path-following algorithm. This computation should be within the capabilities of a modern microprocessor. If such an algorithm could be implemented in the tag reader, we believe it would provide a robust read process and also allow the bit frequencies of the tag to be more closely spaced, possibly even with overlapping hysteresis loops. The reader would deduce the encoded bit pattern from the tag by looking for the presence or absence of a pair of turning points in the path bracketing a bit-encoding frequency.
Our present laboratory implementation of the stabilized sweep is slow because a network analyzer is used as a coherent detector for the null current. There is also some time spent in transferring the data from the network analyzer to the computer implementing the path-following algorithm. However, we do not believe these limitations to be fundamental to the method and an implementation allowing faster stabilized sweep is under study. Observe that with a stabilized sweep, complete information is obtained with only one sweep direction. Table I compares our approach with a recent alternative [22] for reading a passive RFID tag, also based on sensing the response of an oscillator coupled to the tag. Because the two approaches use very different frequency ranges, it is difficult to give a side-by-side comparison. Our approach is probably a bit more expensive to implement in the tag reader, but does give a much more distinct signal for the individual bits encoded by the passive tag. The theoretical analysis presented here has assumed that there is no coupling between individual resonators of the tag itself -only tag-to-oscillator coupling has been treated. In practice, we have observed some coupling between adjacent tag resonators as a second-order effect. Enabling or disabling an encoding bit causes a slight shift in the resonant frequency of adjacent spiral resonators, due to (weak) side-to-side coupling. Again, we do not see this as a fundamental impediment in practice. A passive tag encodes a fixed bit pattern, and a computer optimization could be used to tweak the resonator values to compensate for coupling between adjacent spirals.

VI. CONCLUSIONS
A method for the analysis and synthesis of multiple hysteresis loops in the frequency-tuning curve of an oscillator coupled to a number of external resonators has been presented. The equivalent coupling admittance seen from the coreoscillator terminals has been derived and its particular form of variation has been explained in detail. From the analysis of this function, conditions to avoid the overlapping of hysteresis cycles in the oscillator tuning curve have been derived. A possible application of coupling induced hysteresis in chipless RFID, using an oscillator as a compact reader, has been explored. This application would take advantage of the vertical frequency jumps (associated with the hysteresis phenomenon) to increase the sensitivity to the tag resonances. A proof of concept based on a Colpitts oscillator, with the coupled inductors implemented through planar spirals, has been presented. The coupled system is analyzed extracting the nonlinear admittance function of the core-oscillator, in addition to the coupled admittance, calculated through an electromagnetic simulation. A new numerical method for the experimental characterization of the hysteresis loops, able to pass through their unstable sections, has also been presented, including the mathematical algorithm and the detailed experimental set-up. This should enable an improvement of the shape the hysteresis cycles for an optimum incorporation into the RFID system. Appendix A This section gives a high-level description of an algorithm for tracing a one-dimensional curve.

Zero-Curve
Arc-length Fig. 19. Geometry of path-following algorithm. The zero-curve, starting from x 0 and advancing in arc-length, has crossed from simplex  0 into simplex  1 across the facet common to both simplices. The method to identify such transversal facets is described in the text. Fig. 19 indicates the geometry involved in such a simplicial path-following algorithm. Our previous implementation of a path-following algorithm [15] required a Jacobian matrix of partial derivatives approximated by finite-differencing, which was slow and somewhat noisy. The use here of simplicial methods (which do not require a Jacobian matrix) seems much better suited to experimental work. Starting at the initial point 0 x described above, the algorithm attempts to follow said zero curve until a specific arc-length goal is reached. The user chooses an arc-length goal large enough to include interesting features in the zero curve such as a hysteresis fold-back.
The algorithm constructs a tetrahedron  0 , centered on 0 x with size determined by a user-supplied parameter "grain". For tuning and control voltages in the range [0,10], a typical grain value might be 0.05. A small grain value gives fidelity along the curve, but at some point goes below the resolution of the various DACs. In three-dimensional space a tetrahedron is a simplex -roughly speaking the simplest object that has volume [32], [33]. A simplex in n-dimensional space is determined by n+1 points that do not lie in any sub-space of lower dimension. For example, four co-planer points in R 3 would violate this condition and not qualify as a simplex. The tetrahedron above has four faces determined by any combination of three vertices. In the general case of a simplex, such faces are called facets.
The supporting theory [29] shows that a simplex typically intersects the zero curve in two facets of the simplex or not at all. In other words, moving along the zero curve in the sense of increasing arc-length, the zero curve enters a simplex through one facet and exits it through a different facet. In the supporting literature this is sometimes called the "in door, out door" principle. If the zero curve enters and exits a simplex more than one time, then the simplex is too big. Degenerate cases in which the zero curve intersects an edge or vertex of a simplex are treated in the literature [30] but can largely be ignored in practice because they are so rare.
Solve this 3×3 system and inspect the values of i  to determine if x is interior to the facet. Using this matrix solve, the algorithm inspects each facet of the initial tetrahedron  0 and determines exactly two facets that cross the zero curve (called transversal in the literature).
The algorithm then reports these two facets to the user along with the thee-dimensional coordinates for each point of intersection. In a typical case, one of the intersection point corresponds to a desired direction of travel along the zero curve -for example, to move the tuning voltage in a desired direction. The user specifies the desired direction of travel along the curve by selecting one the two transversal facets of the initial simplex. The algorithm then begins tracing the zero curve. From tetrahedron  0 (in general, simplex  0 ) the algorithm performs a so-called pivot operation to construct the next adjacent tetrahedron ( 1 in Fig. 19) sharing a unique facet. With appropriately chosen tetrahedra (see the discussion in [29], [32]), the pivot step is a very simple computation with vertices. The algorithm writes the coordinates of the exit point of  0 to a file then repeats the computation on  1 to find an exit facet different from the entrance facet. This basic iterative step is then repeated building a "scaffold" of adjacent tetrahedra all hugging the zero curve. Observe that each time a new tetrahedron is adjoined to the scaffold, another measurement is required to find the label of the newly-added vertex. The DACs are stimulated at three voltages given by the coordinates of the new vertex and the VNA is interrogated over a standard measurement interface, such as GPIB to obtain real and imaginary parts of the null current.
The results presented here are for a mapping R 3 into R 2 , but can be generalized. Imagine a reference VCO that is programmed for a waveform more complex than sinusoidalsay containing a fundamental and second harmonic. Such programming would require four quantities -frequency and amplitude of the fundamental, then magnitude of the second harmonic and its phase w.r.t. the fundamental. Recall again that for the autonomous case, the phase reference is determined by the injection signal, hence is not an unknown in the algorithm. Combined with the tuning voltage for the DUT, there would be 5 input quantities. A null detector would then extract real and imaginary parts from the injection current for the fundamental and second harmonic. All four of these numbers would need to be nulled along the zero curve and simplicies along the zero-curve would be labeled with 4 real numbers. The final result, then, would be a mapping R 5 into R 4 . Again, the supporting theory shows that the zero curve of such a mapping is one-dimensional but now moving around in a 5-dimensional space. The algorithm sketched above extends immediately to this case in which tetrahedra are extended to simplices in R 5 .