Systematic Methodology for the Global Stability Analysis of Nonlinear Circuits

A new methodology for the detection of Hopf, flip, and turning-point bifurcations in nonlinear circuits analyzed with harmonic balance (HB) is presented. It enables a systematic determination of bifurcation loci in terms of relevant parameters, such as input power, input frequency, and bias voltages, for instance. It does not rely on the use of continuation techniques and is able to globally provide the entire loci, often containing multivalued sections and/or disconnected curves, in a single simulation. The calculation of Hopf and flip bifurcations is based on the extraction of a small-signal admittance/impedance function from HB and the calculation of its zeros through a geometrical procedure. The method is ideally suited for the investigation of the global stability properties of power amplifiers and other nonlinear circuits. The turning-point locus, associated with either jump phenomena or synchronization, is obtained by taking into account the annihilation/generation of steady-state solutions that is inherent to this type of bifurcation. A technique is also presented for the exhaustive calculation of oscillation modes in multidevice oscillators and oscillators loaded with multiresonance networks. The new methodologies are illustrated through their application to a power amplifier and a multimode oscillator.



Abstract-A new methodology for the detection of Hopf, flip and turning-point bifurcations in nonlinear circuits analyzed with harmonic balance is presented. It enables a systematic determination of bifurcation loci in terms of relevant parameters, such as input power, input frequency and bias voltages, for instance. It does not rely on the use of continuation techniques and is able to globally provide the entire loci, often containing multivalued sections and/or disconnected curves, in a single simulation. The calculation of Hopf and flip bifurcations is based on the extraction of a small-signal admittance/impedance function from harmonic balance and the calculation of its zeroes through a geometrical procedure. The method is ideally suited for the investigation of the global stability properties of power amplifiers and other nonlinear circuits. The turning-point locus, associated with either jump phenomena or synchronization, is obtained taking into account the annihilation/generation of steady-state solutions that is inherent to this type of bifurcation. A technique is also presented for the exhaustive calculation of oscillation modes in multi-device oscillators and oscillators loaded with multiresonance networks. The new methodologies are illustrated through their application to a power amplifier and a multi-mode oscillator.

I. INTRODUCTION
NSTABILITIES are a major concern in power amplifiers and other active circuits [1]- [5], which often exhibit undesired behavior, such as frequency divisions by two, oscillations at incommensurable frequencies or hysteresis [6]- [9]. These phenomena arise through different types of bifurcations [6]- [10], or qualitative changes in the solution stability or the number of steady-state solutions when a parameter is varied continuously [11]- [15]. There are three main methodologies for the prediction of the stability properties at the simulation stage: the Nyquist analysis, applied to the system characteristic determinant [6]- [7], [10], the Nyquist analysis, applied to the normalized determinant function [3], [16]- [17], and pole-zero identification [1]- [2], [18]- [19], applied to a closed-loop transfer function, associated with the circuit linearization about the particular steady-state solution. Only the last two methods are compatible with the use of commercial harmonic balance (HB). On the other hand, the unstable-behavior regions [4]- [5], This work was supported by the Spanish Ministry of Economy and Competitiveness and the European Regional Development Fund (ERDF/FEDER) under the research project TEC2017-88242-C3-1-R. [7], [20]- [24] in terms of relevant parameters, such as the input power, input frequency or bias voltages, can be efficiently determined through bifurcation detection. In comparison with the stability analysis, one advantage of the bifurcation analysis is that the critical poles are purely imaginary or zero at the bifurcation points, so their detection relies on standard simulation methods based on real frequencies.
When using commercial software, the bifurcation analysis can be carried out with auxiliary generators (AGs) operating under certain conditions [4]- [5], [8]- [9], which depend on the kind of bifurcation to be detected. The Hopf and flip bifurcations [20], respectively associated with the onset of incommensurable oscillations and frequency divisions by two, are calculated taking into account that ordinary oscillations are generated/extinguished with amplitude tending to zero [21]. They are detected by introducing a small-signal AG into the circuit, operating at the oscillation frequency, and solving for the parameter values at the limit oscillation condition [20], given by the zero value of the ratio between the current and voltage of the small-signal AG.
The detection of turning points bifurcations, associated with either jump phenomena or synchronization [5], [21], is more demanding. The circuit operates in nonlinear regime with respect to the AG, and the turning-point condition is given by the singularity of the Jacobian matrix of the AG admittance function [21]. In commercial HB, the derivatives in this matrix are obtained applying finite differences to the AG at each iteration of the optimization procedure [5]. This should provide the parameter values for both a zero value of the AG currentto-voltage ratio and a zero value of the determinant of the Jacobian matrix [5], which is numerically involved.
When using commercial software, the bifurcation loci in terms of two significant parameters, 1 and 2 is obtained through optimization plus parameter switching. One of the parameters (1) is swept, solving the particular bifurcation condition through optimization, in terms of the AG variable(s) and the parameter 2, at each sweep step [20], [25]. There can be two different problems associated with this procedure. One is the usual existence of infinite-slope points versus 1, which would demand a manual parameter switching by the user to 2 or the AG frequency or phase, depending on the type of bifurcation. Another problem is the possible existence of disconnected locus curves, which cannot be predicted through continuation techniques (even in in-house software). In this work a new method for bifurcation analysis is proposed, which, unlike previous ones, is systematic. It efficiently provides multiple locus sections, and even disconnected locus curves, with no need for continuation procedures, such as parameter switching. The new method is based on the introduction of a small-signal AG in the HB software. However, this AG is not used to impose a bifurcation condition, solved through optimization, as done in previous works [20], [25], but to extract a small-signal admittance or impedance function. This function is analyzed through contour plots in a dedicated in-house software. This way all the bifurcation points existing for a given value of 1, in terms of the other parameter 2 (and the AG frequency or phase), are obtained simultaneously, regardless of the distance between these points, in terms of the analysis variables. The procedure is extended to obtain autonomous and subharmonic solution curves versus a parameter . The turning-point bifurcations are calculated taking into account that these points are characterized by the annihilation/generation of two steady-state solutions in a solution curve [6], [22].
The new methodology will enable an exhaustive search of free-running oscillations in symmetric multi-device topologies [26]- [29] and circuits loaded with multi-resonance networks, such as reconfigurable oscillators or crystal oscillators [29]. These circuits often exhibit multiple oscillation modes, which are difficult to detect and control, due to their inherent autonomy and the bifurcation phenomena involved in their generation. In previous works [27]- [28], each mode is analyzed in HB individually. Here, the new geometrical method will enable a systematic solution search in given intervals of oscillation amplitude and frequency. Due to their ease of application, we believe that the new methodologies can be very helpful for nonlinear-circuit designers.
The paper is organized as follows. Section II presents the bifurcation detection from dc regime. Section III describes the calculation of multivalued solution curves. Section IV presents the procedures to obtain the flip, Hopf and turning-point bifurcation loci.

II. BIFURCATION LOCI FROM DC REGIME
Due to the usual decoupling of dc signals in microwave circuits, the most common type of bifurcation from dc regime is the primary Hopf bifurcation [6]- [9]. At this bifurcation, and when varying continuously a circuit parameter , a pair of complex-conjugate poles   j crosses the imaginary axis, which gives rise to the onset or extinction of a free-running oscillation.

A. Analysis method
The exhaustive detection of primary Hopf bifurcations is carried out through a methodology in two stages. The first stage is the extraction of a small-signal admittance or impedance function (Ys or Zs) from the microwave-circuit simulation software. This function must be evaluated versus the two analysis parameters, 1 and 2, and the possible oscillation frequency, . The second stage is the simultaneous calculation of the zeroes of this function through a geometrical procedure [30] that does not rely on continuation.
Without loss of generality, the case of an admittance function Ys will be considered. To obtain Ys, a small-signal current source Is will be connected in parallel at a sensitive circuit node, such as a device terminal. The source Is operates at the unknown oscillation frequency . The function Ys is calculated as the ratio Ys = Is/Vs, where Vs is the node voltage at .
To obtain the primary-Hopf locus versus the two parameters 1 and 2, a nested sweep is carried out in 1, 2 and , evaluating the admittance Ys at each step. The  sweep will go from near zero to the transistor maximum operating frequency.
Thus, for each 1, one will have the two-entry function The exported admittance function is split into real and imaginary parts and the data are structured in two matrices: Note that the above equation can have one or more solutions, i.e., there can be one or more Hopf bifurcations for each 1. In a compact manner, the intersection points will be denoted as: One obtains the primary-Hopf locus by sweeping 1 and representing the solutions of (4) in the plane (1,2) or the plane (1,ω). When obtaining multiple intersection points (ω,2) in a certain 1 interval, the locus will be composed by several sections, as in the case of a folded locus, or by two or more disconnected curves. None of this constitutes an analysis difficulty, since the intersection (4) simultaneously provides all the locus points coexisting for 1, regardless of the distances between these points.

B. Application to a power amplifier
The above procedure has been applied to a demonstrator power-combining amplifier, shown in Fig. 1. Note that the purpose of this prototype is not to obtain a very performant design but to exhibit a variety of instability phenomena, enabling an illustration of the capabilities of the new analysis methodologies. The demonstrator is based on the transistor ATF-34143 and is expected at 3 GHz.
The procedure in (1) to (4) has been used to detect the primary Hopf bifurcations of this amplifier. Initially, the auxiliary current source Is, used for the extraction of the admittance function, is connected to the drain terminal of one of the two transistors. The parameters considered are the gatebias voltage VGS and the drain-bias voltage VDS.   Fig. 2(b) shows the primary-Hopf bifurcation locus in the plane (VGS, VDS), with the experimental points at the oscillation boundary superimposed. For VGS below conduction (VGS < -0.81 V), the dc solution is necessarily stable, so the unstable region should correspond to the inside of the locus, as confirmed by the measurements. When connecting the current source Is to the gate terminal one obtains the same Hopf-locus points, superimposed in Fig. 2(b). This is because the bifurcation condition Ys = 0 is fulfilled at all the circuit nodes, so there will be no dependence on the analysis node, provided that the node enables a sufficient observability. One relevant advantage of the new method is that it enables a global exploration of the circuit global stability properties through a simple inspection of the contours Fr(1) and Fi(1). The parameter sweeps should include intervals in which the transistor is necessarily stable, for instance VGS below the conduction threshold or VDS = 0 V. Then, one can ensure that the amplifier is stable for all the possible parameter combinations if no intersection points between the two contours are obtained. Thus, the method is very powerful to assess the small-signal stability of a particular design.
To illustrate this, a stabilization resistor Rs will be connected between the transistor gates. Fig. 2

(d) presents the contours
showing very few intersection points. As the resistor value decreases, the region of these intersection points decreases too and with the standard value Rs = 51 , the contours Fr(1), Fi(1) do not even exist. To illustrate the effect of Rs on the Hopf locus, Fig. 2(b) includes the loci obtained for Rs = 500  (purple) and Rs = 450  (green). The unstable region decreases with Rs. In the experimental characterization Rs = 51  has been chosen to ensure not only stable behavior, but also a sufficient stability margin.
The unstable spectrum obtained experimentally for VGS = -0.6 V, VDS = 1.5 V and Pin = -30 dBm without a stabilization resistor is shown in Fig. 3(a). The stable spectrum obtained experimentally with the stabilization resistor Rs = 51  is shown in Fig. 3(b). With this resistor, and for VGS = -0.7 V, Pin = 3.98 dBm, the amplifier exhibits an output

Lg, wg
Ld, wd power Pout = 10 dBm and an efficiency of 20%. The output spectrum is shown in Fig. 3(c).
Although the instability of the power amplifier could easily be suppressed with the aid of Rs, this resistor will not be used in the subsequent large-signal stability analyses of this amplifier. This is because the aim is to demonstrate the capability of the new methodologies to detect a variety of instability phenomena.

C. Application to a multi-mode oscillator
The method will also be illustrated through its application to the ring oscillator in Fig. 4, susceptible to oscillate at different modes. It is based on the transistor ATF-34143. As will be shown, it exhibits primary Hopf bifurcations around four frequencies: 1.2 GHz, 1.6 GHz, 2.6 GHz and 4.5 GHz.
For the analysis of the primary Hopf bifurcations, a smallsignal current source Is is connected to the drain of one of the transistor devices (Fig. 4). The corresponding loci are calculated versus the gate bias voltage 2 = VGS and the tuning voltage of the varactor diode 1 = VD. Fig. 5(a) presents the intersection points [ωH(η1),η2,H(η1)] obtained through a coarse frequency sweep. The diagram shows three clusters of points, around 1.5 GHz, 2.6 GHz and 4.5 GHz, corresponding to the primary Hopf bifurcations. Thus, the method is capable to predict instabilities at largely spaced frequencies. In fact, it is able to detect all the possible oscillation modes in the interval considered in the frequency sweep. Once the oscillation modes have been detected through the above coarse sweep, the bifurcation loci are calculated through a finer sweep about each detected frequency.
The resulting loci are shown in Fig. 5(b). For each oscillation mode there is a distinct locus curve, so one can immediately know the parameter values (VD, VGS) at which each mode is generated or extinguished. The oscillation frequency corresponding to each mode is indicated, though this autonomous frequency exhibits slight variations through the locus curve. The computational cost is low, since modes at nearby frequencies are obtained simultaneously, such as those at 1.2 GHz and 1.6 GHz. As analyzed later, in this symmetric topology, the phase shift between the transistor stages may be different for the various oscillation modes, though all these modes are detected using a single small-signal current source Is. At low and high VGS, the dc solution is stable for all the VD values. As shown in Fig. 5(b), the primary Hopf loci corresponding to 2.6 GHz and 4.5 GHz exist for all the VD values, and have nearly vertical left and right sections. The primary Hopf locus at 1.6 GHz is composed of two curves. One of the curves exists from VD = 0 V to VD = 2.14 V and the other one exists for VD > 3.86 V. For VD>4.93 V there is also a fourth primary-Hopf bifurcation locus at 1.2 GHz. Each Hopf locus has a left and right section where the corresponding pair of complex-conjugate poles cross to the RHS and to the LHS, respectively, when increasing VGS from low value. The middle regions, comprised between all the direct and inverse Hopf locus sections, contain three pairs of complex-conjugate poles on the RHS for VD<2.14 V, and four pairs of complex-conjugate poles on the RHS for VD>4.96 V.
Due to the points with infinite slope of the locus curves and the fact that the oscillation frequencies detected in Fig. 5(a) are largely spaced, obtaining these loci with the usual bifurcation analysis methods [8], [27]- [28] would be cumbersome. The   [2], [18]. This has been carried out in two stages, initially using a coarse frequency sweep to detect the circuit resonances and then using a fine frequency sweep about each of the detected resonances. The results of these detailed pole-zero identifications have been assembled and are presented in Fig. 5(c) and Fig. 5(d).

GHz
In Fig. 5(c), a constant value VD = 1 V has been considered, increasing VGS from -1.2 V to 0 V. The real part of the dominant poles has been represented versus VGS. When increasing VGS from -1.2 V, the oscillator is initially stable. Then it undergoes a fast sequence of three direct Hopf bifurcations, so, in the interval VGS = -0.86 V to VGS = -0.5 V, the circuit exhibits three pairs of complex-conjugate poles on the RHS. At VGS = 0.5 V the pair of poles at 2.6 GHz crosses to the LHS in an inverse Hopf bifurcation, which is followed by two other inverse Hopf bifurcations, at which the poles at 4.5 GHz and 1.6 GHz cross to the LHS. All the VGS values at which the real part of the poles crosses through zero are accurately predicted by the primary Hopf loci in Fig. 5(b). This can be verified following the arrow corresponding to VD = 1 V. Note that these loci provide the pairs (VGS, VD) at which the oscillator undergoes a primary Hopf bifurcation. In Fig. 5(d), a constant value VD = 9 V has been considered, increasing VGS from ̶ 1.2 V to 0 V. All the direct and inverse Hopf bifurcations are well predicted by the Hopf loci in Fig. 5(b), as can be verified by following the arrow corresponding to VD = 9 V.
As will be shown in the next section, the solutions that are more likely to be observed experimentally are those arising from a stable dc regime or a dc regime with a small number of unstable poles [31]. This is because the unstable poles of the dc solution are transferred to the generated oscillatory periodic solution, in agreement with the well-known bifurcation relationships [22], [32]. Once in the oscillatory regime, and when varying the analysis parameter, these poles must cross to the LHS, in a series of inverse bifurcations, for the solution to become stable. Increasing VGS from low value (following the arrow), the oscillation at 2.6 GHz arises from a stable dc regime. The rest of periodic solutions arise from a dc regime with one or more pairs of complex-conjugate poles on the RHS and are less likely to become stable. On the other hand, for VD < 1.2 V and VD > 3.86 V when decreasing VGS from about 0 V, the oscillation that is generated from a stable dc regime is the one at 1.6 GHz. For VD > 3.86 V, there is a high distance between the Hopf bifurcation at 1.6 GHz and the rest of bifurcations. Thus, one can expect the oscillations at 2.6 GHz and 1.6 GHz to be observed experimentally.
In summary, both Fig. 5(a) and Fig. 5(b) enable an exhaustive exploration of the oscillation modes under variations in two relevant parameters. To our knowledge, bifurcation loci of realistic multi-mode oscillators, with the complexity of those in Fig. 5(b), have never been obtained in any previous work. As shown in the following sections, these loci will be very useful for a reliable design of reconfigurable oscillators.

III. CALCULATION OF SOLUTION CURVES
This section presents a methodology to obtain multivalued and possible multi-section periodic and quasi-periodic solution curves versus an arbitrary parameter , in a systematic manner, without using continuation. The method is an extension of the one described in Section II. However, the admittance function YAG required to obtain the solution curves is nonlinear and must be extracted with a HB simulation, instead of a small-signal analysis. The extraction is carried out with an auxiliary generator [20], which is a voltage generator (playing the role of the oscillation) of amplitude AAG and phase AG, in series with an ideal bandpass filter at its operation frequency AG. The AG is introduced in parallel at a sensitive circuit node [20], such as a device terminal.
In addition to the amplitude AAG, a variable AG is considered, corresponding to the AG frequency, AG = AG, in an autonomous regime, or to the input-source phase, AG = in, in a forced regime (after setting the AG phase to zero, AG = 0). The YAG function is calculated through HB, as the ratio between the AG current and voltage. A triple sweep is performed in the analysis parameter , AG and AAG. When AG = in, the AG sweep is carried out between 0º and 360º. At each AG step, AAG is swept from 0.01 to a few volts (depending on the transistor ratings). A HB simulation is carried out at each sweep step. The inner AAG sweep allows taking advantage of the inherent continuation of the HB inner tier, using the last solution point as initial guess for the next point. One-tone HB is used in the case of a periodic regime and a two-tone HB is used in the case of an autonomous quasi-periodic regime, with a number of harmonic or intermodulation terms as high as required.
All the solution points S() existing for  are exhaustively calculated from: The entire solution curve(s) are given by the whole set of intersection points S() resulting from the  sweep.

A. Application to the multimode oscillator
In case of oscillation modes associated with circuit symmetries (as in the circuit of Fig. 4), the function YAG must be extracted using one AG per transistor stage, in order to sustain the oscillation through the whole structure. The AGs will be connected at equivalent nodes (e.g., at the gate or drain terminals) of the transistor stages. In the presence of circulant passive/active matrixes [34]- [37], the oscillation modes will have the same amplitude at equivalent circuit nodes and different phase shift 2k/N, where k = 1 to N ̶ 1 and N is the number of transistor stages. For each mode, the frequencies and amplitudes of the AGs will be the same [35]- [36], AG and AAG, and the phase shifts between the AGs will be sequentially fixed at 2k/N, where k = 1 to N ̶ 1. Then, the YAG function will be identical at all the transistor stages [34]- [36].
The new method has been applied to calculate the steadystate oscillations of the circuit in Fig. 4. This provides the solution curves in Fig. 6(a) and Fig. 6(c), traced with dots, in terms of VGS, for constant varactor-bias voltage VD = 1 V and VD = 9 V, respectively. The voltage VGS has been purposely used as parameter, instead of VD, to obtain multivalued solution curves, as multivalued sections are commonly encountered about the conduction threshold [38]. In the two cases [ Fig. 6(a) and Fig. 6(c)], the points at which the oscillation curves are originated and extinguished are consistent with the primary-Hopf bifurcations detected in Fig. 5. The distinct solution curves correspond to modes with different inter-stage phase shifts (0º and 120º), indicated in the figure. All the solution curves exhibit turning points or points of infinite slope, which have been obtained without continuation methods.
The accuracy of the systematic method has been validated comparing the curves in Fig. 6(a) and Fig. 6(c) with those obtained through AG optimization/parameter switching [20], superimposed in solid line. Results are overlapped, except in the upper section of the oscillation mode at 1.6 GHz in Fig. 6(a), where the curve exhibits a small loop (expanded view), virtually impossible to obtain through manual parameter switching.
For VD = 1 V, and increasing VGS from low value, the only mode arising from a stable dc solution is the one at 2.6 GHz. This is the only mode that could be observed experimentally, as well as other subharmonic and quasi-periodic solutions resulting from bifurcations of this mode. Measurement points of the mode at 2.6 GHz (before undergoing a frequency division by two, at the point indicated with "") are shown in Fig. 6(a). This mode could be observed experimentally for VGS between -0.945 V and -0.96 V. Fig. 6(b) shows the output spectrum for VGS = -0.95 V, with fundamental output power -6.2 dBm. None of the other oscillation modes detected in Fig. 5 could be observed experimentally, in consistency with the sequence of Hopf bifurcations from dc regime and the discussion at the end of Section II.C. For VD = 9 V, and when decreasing VGS from high value, the only mode arising from a stable dc solution is the one at 1.6 GHz, though in experiment there was a small frequency deviation [ Fig. 6(d)]. We have verified that the spectrum in Fig. 6(d) is not due to a frequency division by two of the solution in Fig. 6(b). Actually, the fundamental frequency in Fig. 6(b) is 2.604 GHz, whereas the fundamental frequency in Fig. 6(d) is 1.4 GHz. Experimental points, corresponding to the 1.6 GHz mode, are superimposed in Fig. 6(c). By varying VD and VGS the circuit behaves as a dual-mode oscillator at 2.6 GHz and 1.4 GHz.

B. Application to a power amplifier
The method in equation (5) has been applied to obtain the oscillatory solution curves of the power amplifier in Fig. 1, versus the input power Pin. When varying Pin, both quasiperiodic solutions, at the input-source frequency in and the autonomous oscillation frequency a, and subharmonic solutions, at in/2, are obtained. To trace the quasiperiodic curves, the variable used in (5) is AG = AG (playing the role of a) and a two-tone HB simulation at in and AG must be performed.
The curves have been traced versus Pin, for VDS = 1.5 V and three different values of VGS. Results obtained for VGS = 0.5 V are traced with dots in Fig. 7(a). For zero input power, the circuit exhibits a free-running oscillation, in consistency with the primary-Hopf locus in Fig. 2(b). When the input power is introduced, the mixing of in with the oscillation frequency a gives rise to an autonomous quasi-periodic regime. Increasing Pin, the quasi-periodic solution is extinguished at the inverse secondary Hopf bifurcation H [6]- [7], [20]. When further increasing Pin, a subharmonic curve arises at the flip bifurcation F [6]- [7], [20]. The subharmonic curve exhibits two turning points: TP1 and TP2. One of them (TP1) is located at a small parameter distance from the quasi-periodic solution curve. No measurement points are superimposed since the represented oscillation amplitude corresponds to the drain voltage of one of the two transistors at a and the odd multiples of the oscillation frequency are ideally cancelled at the output node. Their observation would only be due to discrepancies in the element values of the two amplifier branches. All the solution curves have been validated through AG optimization/parameter switching, with the results superimposed in solid line in Fig. 7(a).  Fig. 7(b) presents the output spectrum for Pin = 0 dBm, corresponding to an autonomous quasi-periodic solution, as predicted by the solution diagram of Fig. 7(a). Fig. 7(c) presents the output spectrum for Pin = 10 dBm, corresponding to a frequency division by 2, in consistency with Fig. 7(a). The analysis in Fig. 7(a) predicts an interval of stable periodic solutions at in between the Hopf bifurcation H and the flip bifurcation F. However, these stable solutions could not be measured for this VGS value. This should be due to the existence of turning points in the experimental quasi-periodic curves, giving rise to an overlapping of these quasi-periodic solutions with the stable periodic solution at in. Those turning points could not be found in simulation for VGS = -0.5 V, which is attributed to modelling inaccuracies. However, when changing the transistor gate-drain Cgd capacitance, as a small dispersion test, from 0.15 pF (nominal value) to 0.25 pF, turning points arise in the quasiperiodic curve [ Fig. 7(d)], which would lead to jumps from quasiperiodic to subharmonic regime and vice versa.
For VGS = 0.37 V and VGS = 0.6 V one obtains the curves in Fig. 8(a), traced, respectively, with blue and red dots. For VGS = 0.37 V, the turning point TP2 of the subharmonic curve has vanished and this curve exhibits a loop before decaying to zero at the flip bifurcation F, as in the case of Fig. 7(d). In this case, there is a Pin interval for which the subharmonic curve does not coexist with the stable periodic solution at in (comprised between H and F). This situation, which is not found for VGS = 0.5 V and VGS = 0.6 V, should facilitate the experimental observation of this periodic solution. See the output spectrum obtained for Pin = 18 dBm in Fig. 8(b). For VGS = 0.6 V, the subharmonic solution is composed by two disconnected curves, which would be difficult to detect with an ordinary method. In view of this result, one can suspect that bifurcation diagrams obtained through continuation will often miss curve sections. Note that there is no way to predict the frequency division in the closed curve through a stability analysis of the periodic regime at in. After its detection with the new methodology, this solution curve has also been obtained through AG optimization and parameter switching (in solid line).
It can be of interest to compare the advantages and drawbacks of the new method (M1) and the one based on AG optimization/parameter switching (M2). When tracing a nonproblematic curve section, the computation time of M2 is shorter, since it directly addresses each solution point, whereas M1 is based on the extraction of the nonlinear function ( , ) AG AG AG YA   . However, M2 demands user surveillance and manual parameter switching, so the total time devoted to the circuit simulation may be very long in complex solution curves. Furthermore, there can be situations in which no convergence is obtained with M2. This is because the two-tier representation (where YAG = 0 is the outer tier and the pure HB system the inner tier) limits the parameter switching to only three variables, though the nonlinear circuit is actually governed by a HB system with many more state variables. These convergence failures leave the user without any possibility to complete the solution curves. In addition, M2 is not exhaustive, unlike M1, so it will miss isolated curves, like the one in Fig. 8(a).

IV. BIFURCATION LOCI FROM PERIODIC REGIME
The three types of local bifurcations from periodic regime are flip, Hopf and turning-point bifurcations [6]- [10]. At flip and Hopf bifurcations an oscillation is generated/extinguished with zero amplitude. At turning-point bifurcations, the solution curve exhibits an infinite slope. The exhaustive calculation of the corresponding bifurcation loci in a two-parameter space (1,2) is presented in the following paragraphs.

A. Hopf and flip bifurcations
At a Hopf bifurcation from a periodic regime at in, a pair of complex-conjugate poles σ ± jωa, where ωa is incommensurable with in, crosses the imaginary axis of the complex plane, giving rise to the generation or extinction of a quasi-periodic regime [6]- [10]. On the other hand, at a flip bifurcation from a periodic solution at in [6]- [10], a pair of complex-conjugate poles σ ± jωin/2 crosses the imaginary axis of the complex plane, giving rise to the generation/extinction of a subharmonic oscillation at ωin/2. The Hopf (flip) bifurcations can be detected introducing a small-signal AG at ωa (in/2) into the circuit. Note that the circuit operates in large-signal conditions with respect to the input source at in but in small signal with respect to the AG, since its amplitude is set to a small value. This AG will be used to obtain the small-signal admittance YAG, which should be zero at both the Hopf and flip bifurcations.
For each 1, a double sweep is performed in AG and 2, where AG = ωa, for the Hopf-bifurcation detection, and AG = in, for the flip-bifurcation detection. The double sweep in AG and 2 has a very low computational cost due to the small-signal amplitude of the AG. All the bifurcation points B(η1) existing for 1 are exhaustively calculated from: The entire bifurcation curve(s) are given by the whole set of intersection points B(η1) resulting from the 1 sweep.

B. Turning points
At a turning point bifurcation, obtained when varying a parameter , two points of a same solution curve (S1 and S2) merge and annihilate, which gives rise to an infinite slope. The bifurcation relationship, in a ball of the state-variable space, with radius tending to zero, can be expressed as [6], [22]: where the superscript indicates the number of unstable canonical poles and the symbol  indicates no solutions.
Since a real pole passes through zero at the turning point bifurcation, the difference between the number of unstable poles of S2 and S1 is necessarily n = 1. The relationship (7) is fulfilled at turning points in solution curves of any kind of regime, i.e. dc, periodic or quasi-periodic.
When using the new geometrical method to obtain solution curves, the turning points can easily be detected in the following manner. One should keep track of the number of steady-state solutions obtained for each parameter value (k), expressed as m[η(k)], where k is a counter of the parameter values for which solutions are found. Taking (7) into account, turning points are obtained at the  values where m increases or decreases in 2.
Thus, at a turning-point, the following condition is fulfilled: The above condition enables a straightforward calculation of the turning-point locus, as shown later in this section. Additionally, at the design stage, one can predict the occurrence of turning points by evaluating the Jacobian matrix of the outertier admittance function YAG with respect to AAG and AG. As demonstrated in [21], this Jacobian matrix becomes singular at turning points, which fulfill: By continuity, the function det(YAG) will have opposite sign at the two points that arise or annihilate at the turning point (just before or after this bifurcation). The determinant det(YAG) is easily calculated through a spline interpolation of The ±1 is due to the fact that the determinant has opposite sign at the two steady-state solutions that appear/disappear at the turning point. Thus, both the number of solutions with det(YAG)>0 (m P ) and the number of solutions with det(YAG)<0 (m N ) must increase or decrease in one.
The new procedure for the calculation of turning points will be illustrated through its application to the closed subharmonicsolution curve in Fig. 8(a)  and det η (ϕin,AAG) = 0, which in the case of Fig. 8(a) and Fig. 9 occurs only for two ϕin,AAG values: ϕin = 179º and AAG = 1.504 V, ϕin = 192º and AAG = 1.49 V. Note that obtaining the contour det η (ϕin,AAG) = 0 is not necessary for the evaluation of the turning-point conditions (8) and (10). However, its calculation has very low computational cost and can enable an early detection and suppression of turning points during the design procedure, in a similar manner to what was done in Fig. 2(d) with the stabilization resistor. Conditions (8) are easily extended to an analysis in terms of two parameters 1, 2. This is done by performing two nested sweeps, in 1, 2. At each step of the outer sweep in 1, one stores all the 2 points, denoted as 2T, that fulfill (8). The turning point locus is obtained by representing 2T versus 1.

C. Application to a power amplifier
The bifurcation-analysis method has been applied to the power amplifier in Fig. 1. The bifurcation loci have been calculated in terms of 1 = VGS and 2 = Pin for two values of the drain-bias voltage: VDS = 1.5 V, in Fig. 10(a), and VDS = 3 V, in Fig. 10(b). These two values have been selected since, as shown later, they give rise to qualitatively different bifurcation loci in the plane (VGS, Pin).
For VDS = 1.5 V one obtains a flip locus composed of two disconnected curves (Flip Locus 1 and Flip Locus 2) that cannot be easily detected using ordinary techniques. When VGS is swept from -0.7 V to -0.3 V, the new method simultaneously provides points in the two flip-locus curves. This is shown in Fig. 12(a), which presents the intersections between 1 ,2 ( , ) 0 The versatility of the new method enables an investigation of the effect of an additional (third) parameter on the locus curves. In Fig. 12(b), VDS has been added as a third parameter. Increasing VDS, the two disconnected flip locus curves approach each other and merge into a single one. To our knowledge, no analysis of such complexity has ever been performed in combination with HB commercial software. At VDS = 1.5 V and VDS = 3 V, one respectively obtains two flip loci and one flip locus, leading to a qualitatively different behavior. This is why these VDS values have been chosen in this investigation.
As shown in Fig. 10(a), the Flip Locus 1 exhibits a pronounced minimum at VGS = -0.88 V because for that VGS value the frequency of the (stable) dominant complex-conjugate poles approaches the input frequency divided by two. The Hopf locus only exists on the right side of the flip locus. The Hopflocus point corresponding to Pin =  dBm (in the absence of input power) agrees with the primary-Hopf bifurcation obtained in the analysis of Fig. 2 for VGS = -0.83 V and VDS = 1.5 V. The Hopf and flip loci merge at the two co-dimension two bifurcations [12]- [13] indicated in the figure as CT1 and CT2.
The periodic solution at the input-drive frequency in is unstable inside the flip locus and inside the Hopf locus.
To complete the stability portrait, the locus of turning points in subharmonic regime has also been represented in the plane (VGS, Pin) [see Fig. 10(a)]. Performing a double sweep in VGS and Pin, one obtains two turning-point locus sections, shown in Fig. 10(a), respectively denoted as TP1 Locus and TP2 Locus. In fact, the turning points observed in the subharmonic solution curves of Fig. 7 and Fig. 8 (also denoted as TP1 and TP2) belong to these two sections. Closed subharmonic solution curves [such as the one in Fig. 8(a)] are obtained in the region of the plane in Fig. 10(a) located between the two sections of the flip locus. After the intersection of the TP2 Locus and the Flip Locus 2 one obtains open subharmonic curves, as the one in Fig. 7(a).
Turning-point bifurcations may correspond to two different phenomena: jumps and synchronization [20]- [21]. When departing from a quasi-periodic solution and varying any of the two analysis parameters in its full range, synchronization points will be those traversed without ever crossing the Hopf locus. This is because synchronization is the only mechanism for the quasi-periodic solution to transform into a subharmonic one. When also traversing the Hopf locus, synchronization points may still exist but their identification is demanding, since it requires the detection of the global saddle-connection bifurcation [12], [20], [31]. The size of the parameter intervals for which this bifurcation occurs is typically very small and located near the co-dimension two bifurcations. Thus, this bifurcation is not generally very relevant.  From the previous discussion, one can assure that turning points obtained for Pin values between CT1 and CT2 will be synchronization points. However this interval will actually be broader. Due to the occurrence of this synchronization phenomenon, subharmonic solutions exist for Pin values smaller than the one at which the flip bifurcation takes place.
Measurement points at the boundary between the parameter regions with stable behavior at in and subharmonic behavior at in/2 are superimposed with squares in Fig. 10(a). For VGS values on the right side of the Hopf locus, as Pin increases from low value, one initially observes a quasiperiodic spectrum and then a frequency-division by two. The experimental points where the transitions occur are represented with diamonds. At these transitions, a triangular spectrum, characteristic of synchronization, is observed (Fig. 11), which is in agreement with the presence of the turning-point locus of Fig. 10(a). As shown in Fig. 10 there is good qualitative agreement between simulation and measurements. Discrepancies are due to the different shape of the experimental loci, which is attributed modelling inaccuracies. In the next subsection, the bifurcation loci will be validated through pole-zero identification, which is an independent analysis method, relying on identical models of the circuit components. Thus, there should be a perfect correspondence between the two methods.

D. Validation with pole-zero identification
The predictions by the flip and Hopf loci have been validated through pole-zero identification. In Fig. 13, the real part of dominant poles obtained for VDS = 3 V and two input-power values, Pin = 4 dBm and Pin = 17.95 dBm, has been represented versus VGS. For Pin = 4 dBm [ Fig. 13(a)] and low VGS, the circuit is initially stable, as expected with the transistor below conduction. There is a stable pair of dominant complexconjugate poles σ ± jωa, where ωa/(2)  1.4 GHz. When varying the parameter, ωa will vary continuously. Increasing VGS, the pair of poles at σ ± jωa transforms into two distinct pairs of complex-conjugate poles at in/2 at VGS = 0.98 V.
These two distinct pairs of poles are denoted as σ ± jωin/2 and σ′ ± jωin/2. The pair of poles σ ± jωin/2 crosses to the RHS at VGS = 0.94 V, in a direct flip bifurcation (F1), from which the periodic solution at in becomes unstable. This flip bifurcation is consistent with the prediction by the flip locus in Fig. 10(b), as verified by following the arrow at Pin = 4 dBm, indicating the VGS variation. Every crossing of the arrow through the flip locus corresponds to a flip bifurcation, at which the pair of poles σ ± jωin/2 crosses the imaginary axis. At VGS = 0.79 V the same pair of poles σ ± jωin/2 crosses to the LHS in an inverse flip bifurcation (F2) that is also accurately predicted by the locus in Fig. 10(b). Just after crossing (VGS = 0.77 V), the two pairs of unstable complex-conjugate poles σ ± jωin/2 and σ′ ± jωin/2 merge into a pair of complex-conjugate poles σ ± jωa and cross to the RHS in a Hopf bifurcation (H). This Hopf bifurcation is also consistent with the prediction by the Hopf locus in Fig. 10(b), as verified by further following the arrow corresponding to Pin = 4 dBm, as VGS increases. The crossing of the arrow through the Hopf locus indicates a Hopf bifurcation, at which the pair of poles σ ± jωa crosses the imaginary axis. There is only a small VGS interval with stable behavior, in consistency with the neighborhood to the co-dimension two bifurcation CT1' predicted by the loci in Fig. 10(b). For Pin = 17.95 dBm [ Fig. 13(b)] and low VGS, the circuit is initially stable. There is a flip bifurcation at VGS = 0.88 V (F1), where a pair of complex-conjugate poles σ ± jωin/2 crosses to the RHS. At VGS = -0.45 V , this same pair of poles crosses the imaginary axis to the LHS, in an inverse flip bifurcation (F2) that is also consistent with the prediction of the locus in Fig. 10(b), when following the arrow corresponding to Pin = 17.95 dBm. Just after crossing (VGS = -0.445 V), the stable pair of poles σ ± jωin/2 merges the second stable pair σ′ ± jωin/2 and they become a pair of complex-conjugate poles at an incommensurable frequency, σ ± jωa, which crosses to the RHS at VGS = 0.39 V in a direct Hopf bifurcation (H). Note the consistency between the results of the two methods (by following the arrow corresponding to Pin = 17.95 dBm), even though the bifurcations are closely spaced. Identical consistency has been found in all the additional validations performed by varying the input power and biasing conditions.

V. CONCLUSION
A conceptually simple methodology for an exhaustive bifurcation analysis in nonlinear circuits, which does not rely on continuation techniques, has been presented. Its main advantages are its systematic quality and its global exploration capabilities. When tracing the bifurcation loci, the method simultaneously provides bifurcation points corresponding to different locus sections and even disconnected locus curves. The method is intended for application in combination with commercial harmonic-balance software, from which an admittance function, depending on the analysis parameters, is extracted. It enables an efficient calculation of instability boundaries in nonlinear circuits prone to instability, such as power amplifiers. It also allows a systematic investigation of the complex problem of multiple oscillation modes in multidevice oscillators and oscillators loaded with multi-resonant networks. The calculation of the turning point locus has a high computational complexity with ordinary methods, which impose a singularity condition to the nonlinear-system Jacobian matrix. Here it is obtained in an original manner, taking into account the annihilation of solution points occurring at this type of bifurcation. The new method has been successfully validated through its application to a power amplifier and a multi-mode symmetric oscillator.