Analytical and Numerical Bifurcation Analysis of Circuits Based on Nonlinear Resonators

—This work presents new frequency-domain methods for the analysis and simulation of circuits based on a nonlinear resonator, with operation ranges delimited by turning points. Insightful analytical conditions fulﬁlled at these turning points are derived, which will enable an identiﬁcation of the effect of each circuit element on their location in the solution curve. The cusp points, or co-dimension two bifurcations at which two turning points merge into one, thus delimiting the multivalued intervals, are directly calculated for the ﬁrst time to our knowledge. In addition, a new numerical method, compatible with the use of commercial harmonic balance, is presented for the straight-forward tracing of the multivalued solution curves, together with a new procedure to determine the locus of turning points in terms of any two analysis parameters. This relies on the use of a new mathematical condition to obtain a surface of turning points in the space deﬁned by the two parameters and the input power. The methods have been applied to a wireless power-transfer system based on a recently proposed conﬁguration, obtaining very good agreement with the experimental results.


I. INTRODUCTION
R ECENTLY, interesting applications of nonlinear res- onators to improve the behavior of energy harvesters and wireless power-transfer systems have been proposed [1]- [7].In a nonlinear resonator, the oscillation frequency depends on the amplitude [8], [9], and when excited with a forcing signal, there is a folding of the resonance curve that typically exhibits two turning points (TPs), or infinite-slope points.As a result of this folding, the response exhibits a broader frequency bandwidth (BW) in comparison with that of a linear resonator (see Fig. 1(a), reproduced from [6]).The works [1] and [6] propose the application of this BW enhancement in RF energy harvesters.This is because the high voltage levels required for a proper rectification can be achieved with improved tolerance to the frequency detuning [6] caused by aging, coupling effects, etc.In the case of Fig. 1(a) [1], [6], Fig. 1.Multivalued response of a nonlinear resonator [6].The turning points at which physical jumps are obtained when varying input signal are indicated.(a) Variation of the output amplitude versus the input frequency.(b) Variation of the output amplitude versus the input amplitude.the higher-amplitude turning point (TP2) [10] determines the maximum operation frequency, since when reaching this point (with an infinite slope) a (discontinuous) physical jump takes place to the lower-amplitude section of the multivalued curve.
When varying the excitation amplitude at a constant frequency, one obtains a significantly flat higher-amplitude section [Fig.1(b)].As demonstrated in [5] and [7], this is useful in near-field resonant-based wireless power-transfer systems [11], [12].In those systems, the power entering the second resonator is highly dependent on the coupling factor, which, in turn, varies with the distance between the transmission and reception coils and/or misalignments.As demonstrated in [3]- [7], the sensitivity to the coil positions can be reduced by taking advantage of the flattened response of the nonlinear resonator [Fig.1(b)], delimited by the turning point (TP2).
The objective of this work is the development of new frequency-domain methods to facilitate the design and simulation of circuits based on a nonlinear resonator.Attention will be paid to the turning points [8], [10], [13] (corresponding to D-type bifurcations [14], [15]), as they delimit their usable ranges.We will consider two illustrative cases: a single nonlinear resonator, as in the RF energy harvester of [1] and [6], and a linear resonator coupled to a nonlinear one, as in the wireless power-transfer system of [7].Insightful analytical conditions fulfilled at the turning points will be derived, which will enable an identification of the effect of each circuit element on their location in the solution curve.The turning point condition can be used for a direct prediction of the multivalued interval in terms of an analysis parameter, such as input frequency, power, or distance.On the other hand, the cusp points [8], [10], [13], or co-dimension two bifurcations at which the two turning points of the solution curve merge into one, provide the limit of the multivalued interval.Here these cusp points will be directly calculated for the first time to our knowledge.
Besides the analytical investigation, a novel numerical methodology, compatible with the use of commercial harmonic balance (HB), will be proposed for the straightforward tracing of the solution curves and turning point loci.Note that an ordinary parameter sweep is unable to pass through the turning points, as these are singular points of the HB system [14]- [17].There can be a loss of convergence in its neighborhood or a discontinuous jump to another section of the multivalued curve, so the empirical estimation of the turning points based on these jumps is unreliable, in general.In commercial HB, one way to pass through the turning points is to introduce an auxiliary generator (optimized to fulfill a nonperturbation condition) and change manually the sweep parameter when convergence fails, as the turning point is approached [16]- [19].However, this manual parameter switch becomes demanding when numerous tests are required.For a more efficient analysis, a contour-tracing method has been recently proposed [20], [21], which involves a two-step procedure.The first step is the scattering-parameter simulation of the input network to calculate its Norton equivalent at the fundamental frequency.The second step is the calculation of a nonlinear current function with the aid of an auxiliary generator [16], [17].To trace the solution curve, one must combine the results of the two independent analyses, which can be demanding and prevent the quick feedback required in the design/optimization process.
Here, a new analysis method is presented, which allows obtaining the solution curves and turning points loci in a single simulation with the aid of some equations easily introduced in commercial HB.The method is extended to obtain, for the first time to our knowledge, the locus of turning points versus any pair of arbitrary parameters without optimization/continuation procedures.This relies on the use of a new mathematical condition to obtain a surface of turning points in the space defined by the two analysis parameters η 1 , η 2 and the input power P in .Then, the intersection of this surface with planes defined by distinct P in values will provide the corresponding turning point loci in terms of η 1 and η 2 .This turning point surface enables an exhaustive calculation of the family of turning point loci resulting from the variation of the input power.The new analytical and numerical methods are expected to be useful to researchers on novel interesting applications of circuits based on a nonlinear resonators, such as those in [1]- [7].
This article is organized as follows.Section II presents the analytical investigation of a single nonlinear resonator with emphasis on the turning point mechanism.Section III presents an analogous investigation in the case of a linear resonator coupled to a nonlinear one.Finally, Section IV presents the new numerical methods to obtain the multivalued solution curves, as well as the locus of turning points in terms of any two analysis parameters.

A. Turning Point Condition
Let the single nonlinear resonator in Fig. 2 be considered.The nonlinear capacitor is represented with the describing function Q(V ) associated with its instantaneous charge q(v), where v(t) is the excitation voltage.Applying Kirchoff's laws at the excitation frequency ω, one obtains where φ is the opposite of the phase shift (acting as a state variable) between the node voltage and the input current I g , and V is the voltage amplitude.On the other hand, H (V, ω) is a nonlinear function that compactly describes the circuit response to the voltage amplitude V at the frequency ω at the analysis node.In the presence of the input source I g at ω, the circuit solutions must fulfill H (V, ω) − I g e j φ = 0.This complex equation is composed of two real equations in the two variables V and φ.At the turning points, the Jacobian matrix associated with (1) becomes singular, which is also the condition for a bifurcation of D type [14]- [17].This Jacobian matrix is calculated with respect to the two state variables of (1), given by V and φ where the superscripts r and i indicate real and imaginary parts.One can also solve for I g sin φ and −I g cos φ using (1), which provides det Note that the above condition only depends on the excitation amplitude V and frequency ω.Using the definition of H given in (1), the above determinant can be written as To obtain the locus of turning points in the plane defined by ω and V , one should solve the scalar equation ( 4) in terms of the voltage amplitude V for each ω, which will provide the turning points denoted as V T .Using (1), the input-current amplitude at each turning point is simply calculated from (5) To distinguish the nonlinear effects from the linear ones, the charge describing function Q(V ) will be expressed in terms of a linear term CV and a purely nonlinear function Replacing ( 6) into (4) and separating the terms associated with the linear resonance of C and L (under the conductance G) from the nonlinear ones, one obtains the singularity condition Note that V is an amplitude, so V > 0. The above turning point condition is composed of three terms.The first term, purely linear, can only be positive, so the zero-value condition requires a negative contribution.If the nonlinear component of the charge function Q NL (V ) and its voltage derivative ∂ Q NL /∂ V have the same sign for all V , the third term will also be positive.In these conditions, only the second term can be negative, and the folding to either higher or lower frequencies than the linear resonance frequency ω o = 1/ √ LC will depend on the sign of Q NL (V ) and ∂ Q NL /∂ V .If their sign is positive (negative) for all V , all the nonlinear-resonance curves will fold leftwards (rightwards), that is, toward lower (higher) ω values, which will be due to the dominant contribution of −1/(Lω) (the dominant contribution of Cω).From the inspection of (7), both for a larger resistive effect (higher G) and a larger detuning with respect to ω o , one will have a larger positive first term and, thus, will need a stronger negative effect of the nonlinearity to get det(V, ω) = 0.

B. Polynomial Approximation
As shown in [3]- [7], when limiting the Taylor series expansion of the junction capacitance to the second order, the charge of two antiseries diodes can be described as q(t) = C o v(t)/2+ q 3 v 3 (t), where C o is the zero-bias junction capacitance of the individual diodes and q 3 may be positive or negative, depending on the diode characteristics.The associated describing function is given by For q 3 < 0, one has Q NL (V ) < 0 and ∂ Q NL /∂ V < 0 for all V.In these conditions, and taking (7) into account, the resonance curves will fold rightwards (to ω > ω o ), as this provides a negative value of the second term.This result, in agreement with [4], has been validated with a numerical example based on the diode SMV1470-004LF.The back-to-back connection of two of these diodes has been approximately modeled with the coefficients C = 55 pF and Q 3 = −2 × 10 −12 C/V 3 .The rest of element values are shown in the caption of Fig. 2. The solution curves calculated through (1) for different input-current amplitudes, in terms of V versus ω, are shown in Fig. 3(a).As expected, all the curves fold rightwards.For the highest input current, the analytical results have been compared with default HB simulations performed by sweeping the input frequency ω.Because the aim is to validate the analytical results (based on the describing function), only the fundamental frequency is considered.Because these default HB simulations are unable to circumvent the turning points, a discontinuous jump to the lower section of the curve is obtained at the first turning point.The rest of the two curves are overlapped.Note that Fig. 3(a) also includes the curves obtained with the HB-based method in Section IV, able to pass through the two turning points.Instead of directly addressing the circuit with its input current source at a particular value I go , this method is based on the calculation of the nonlinear function H (V, ω) with the aid of an auxiliary generator introduced in commercial HB.Once the function H (V, ω) is available, the solution curve for any I g is given by the corresponding contour level of the surface |H (V, ω)|.Thus, we can obtain a whole family of solution curves (for distinct I g values) from the same function H (V, ω).The contours can have any shape, so they naturally provide multivalued curves versus ω, with no need for continuation methods such as parameter switching [17], [20].Section IV provides the necessary implementation details.The resulting numerical curves are fully overlapped.
Fig. 3(b) presents the resonance curves obtained for R = 1/G = 1 k, at different input-current amplitudes.In agreement with the analytical derivations, the curve folding is observed from a frequency closer to the original resonance one ω o = 1/ √ LC of 12 MHz, and the required input amplitudes are smaller than those in the curve family of Fig. 3(a).For the highest input current, the analytical results have been compared with default HB simulations.The curve family has also been compared with that obtained with HB-based method in Section IV and the results are overlapped.
When considering the charge function Q(V ) = CV +Q 3 V 3 , one can explicitly obtain the locus of turning points by replacing (8) into (4) and solving for V in terms of ω.This leads to the following expression for the voltage amplitude at the turning points: From the inspection of (9), there are two turning points at each ω, except at the cusp point where the radicand becomes zero.The turning point locus given by (9) has been superimposed on the solution curves of Fig. 3(a) and (b) in solid black line.This locus is the curve that passes Fig. 3. Analysis of the nonlinear resonator in Fig. 2. (a) Nonlinear resonance curves for different values of the input-current amplitude, in terms of V versus ω, for G = 1/500 = 2 mS.Default HB simulations are superimposed for I g = 9 mA.HB results obtained with the numerical method in Section IV are also presented for comparison.The locus of turning points obtained through ( 9) is superimposed and the cusp resulting from ( 11) and ( 12 through all the infinite-slope points of the nonlinear-resonance curves obtained when varying the input-current amplitude I g .As stated, the locus is given by (9).The sign of the first term in the numerator must be negative for Q 3 < 0 since only positive values of the expression in (9) are valid, so the turning points necessarily occur for ω > ω o .The negative sign of the square root in (9) provides the turning points in the upper section of the nonlinear-resonance curves, which, when increasing the input frequency, give rise to a downward jump to the lower section.The positive sign of the square root provides the turning points in the lower section, which, when reducing the input frequency, give rise to an upward jump to the upper section.When the input current I g decreases and approaches the value for which the cusp point is obtained, the downward and upward jumps occur at progressively closer frequency values, which degenerate to the same value at the cusp point.In the curves obtained for lower input currents than the one corresponding to the cusp, there are no turning points, and thus no physical jumps.Note that in the case Q 3 > 0, the sign of the first term in the numerator of ( 9) must be positive, so the turning points will necessarily occur for ω < ω o .The positive (negative) sign of the square root will provide the turning points in the upper (lower) section of the nonlinear-resonance curves.
As gathered from ( 9), the amplitudes V T scale down with √ |Q 3 |, that is, they decrease for a more nonlinear Q(V ).For a higher G (lower quality factor), the turning points will start from a larger frequency detuning, required for a positive radicand.Note that each turning point, given by ω and V T , is obtained for a different input current, which will be calculated replacing ( 9) into ( 5).This provides Thus, one can easily obtain the locus of turning points in the plane defined by ω, I g,T .
As stated, the radicand of ( 9) becomes zero at the cusp point, which is obtained for the frequency In this limit case, the upper and lower turning points of the solution curve merge into one, corresponding to a cusp.This is a co-dimension two bifurcation [8], [10], requiring the fine tuning of two parameters, so, in practice, the circuit can only operate in the neighborhood of the cusp.This cusp has a relevant meaning since for Q 3 < 0 (Q 3 > 0) it provides the minimum (maximum) frequency at which the resonance curves exhibit folding [see Fig. 3(a)].For ω < ω c (ω > ω c ), there are no turning points.For ω > ω c (ω < ω c ), there will be two turning points, each for a different V T (and thus, a different input current).As gathered from (11), the frequency of the cusp is independent of the nonlinearity.On the other hand, the detuning of ω c with respect to ω o will be larger for a higher ratio G/C.
The voltage amplitude at the cusp points is This voltage will be higher for a larger detuning of ω c with respect to ω o .The cusp resulting from (11) and ( 12) is indicated in Fig. 3(a) and (b).Finally, Fig. 3(c) presents the evolution of the turning points loci when varying the conductance G at constant C. The evolution of cusp when sweeping G, taken as an implicit parameter, is traced with a dashed line.

III. LINEAR RESONATOR COUPLED TO A NONLINEAR ONE
The next analytical investigation tackles the case of a series linear resonator coupled to a nonlinear parallel one, as proposed in [7] for wireless power transfer.The schematic is shown in Fig. 4. Note that a low-power system has been implemented here since the focus of this work is the demonstration of the new analysis and simulation methods, instead of the system performance.

A. Equation System
Applying Kirchoff's laws at the excitation frequency ω, one obtains the following set of four complex equations in the four complex unknowns V 1 , V 2 , I 1 , and I 2 : where E g is the input-voltage amplitude, R s is the source resistor, L p and L s are the primary and secondary inductors, respectively, M is the mutual inductance is the nonlinear charge, and G L = 1/R L is the output conductance.The first three equations ( 13)-( 15) are linear with the unknowns I 1 , V 1 , I 2 , V 2 , and are easily solved for I 2 in terms of V 2 and E g .Replacing the resulting expression for 16) one obtains a complex equation depending on V 2 only, which can be expressed on terms of the magnitude |V 2 | and the opposite of the phase shift between V 2 and E g , given by φ.This final complex equation is the following: where φ is the opposite of the voltage phase shift with respect to the input source, and A 1 (ω) and A 2 (ω) are passive linear coefficients depending on ω only, whose explicit expressions are shown in (18).Both |V 2 | and φ are the unknowns of the complex equation (17), which, in an explicit manner, is given by The inspection of the above expression evidences the combined action of two resonators.For k = 0, the secondary resonator becomes isolated, with the only solution |V 2 | = 0 and the primary becomes an ordinary series resonator at ω p = 1/(L p C p ) 1/2 .For 0 < k < 1, the signal from the input source E g reaches the nonlinear parallel resonator through the coupling effects.

B. Frequency Response in Linear Conditions
To get some insight into the frequency response of the circuit in Fig. 4, we will initially replace the nonlinear charge The results of this linear analysis will serve as a comparison reference for the nonlinear case, treated in Section III-C.After replacing (18), the linear transfer function from E g to V 2 , in terms of the complex frequency s, is given by where the denominator is the characteristic polynomial [22] of the coupled circuit.In this polynomial the terms depending on k have been separated from the rest, which is expressed explicitly as a product of the characteristic polynomials of the two individual resonators in the absence of coupling (k = 0).We will now analyze the evolution versus the coupling factor k of the poles of H L (s), given by the roots of the denominator A relevant case is the one in which L p = L s = L and C p = C, and, prior to the coupling, the two resonators have the same resonance frequency ω o = 1/ √ LC.For sufficiently low G L and R s at k = 0, the circuit exhibits two pairs of complex-conjugate poles with near identical imaginary parts and different real parts depending on G L and R s , respectively.This is shown in Fig. 5(a), corresponding to L p = L s = 3.2 μH and C p = C = 55 pF.When k increases, the two pairs of complex-conjugate poles approach each other and, at a given coupling factor k c , they fold in opposite directions.From k c the frequency originally associated with the primary (secondary) resonator decreases (increases).
The results of the pole analysis are consistent with the variation of the magnitude of H L ( j ω) represented versus frequency in Fig. 5(b).The results obtained with the analytical expression (19) are overlapped with those resulting from a numerical linear analysis in commercial HB.In agreement with [11] and [12], two maxima are obtained as k increases.As predicted by the pole analysis in Fig. 5(a), the upper (lower) resonance shifts rightwards (leftwards), and the upper resonance shift is more pronounced than the lower one.Because the higher-frequency poles move away from the imaginary axis, the resonance becomes less pronounced, as seen in Fig. 5(b).On the other hand, the lower frequency poles slightly approach the imaginary axis, and the resonance keeps pronounced.Note that for k = 1, the characteristic polynomial becomes

C. Nonlinear Operation
The complex equation (17), describing the linear resonator coupled to a nonlinear one, can be expressed in terms of a compact nonlinear function H (|V 2 |, ω) by dividing both terms of this equation by A 1 (ω).This provides This way the circuit is formally described in the same way as the single nonlinear resonator of Section II [see (1)].For a given input amplitude E go , the nonlinear resonance curve versus the frequency ω is calculated obtaining the absolute value of H (|V 2 |, ω) and, thus, making the phase shift φ disappear At each ω, the real equation ( 23) provides one or more solutions in terms of |V 2 |.Note that for each pair of values ω, |V 2 | there will be a different phase shift φ that can be calculated replacing these values in (22).
A very simple way to get a whole family of resonance curves for as many E g values as desired is to perform a double sweep in |V 2 |, ω and obtain the surface |H (|V 2 |, ω)|.Then, the nonlinear resonance curve for any E g = E go is given by the contour level This is the procedure followed in the HB method described in Section IV, though in the HB calculation one can consider any number NH of harmonic terms.   in Fig. 6(c) for C p = 41.5 pF, considering the coupling factor k = 0.15 and three different excitation frequencies 15, 15.5, and 16 MHz.The analytical results are compared with those obtained with the HB method of Section IV and with default HB simulations.Because the default HB simulations are unable to circumvent the turning points, a discontinuous jump is obtained at the first turning point.The rest of the default-HB curve is overlapped in all cases.Fig. 6(d) presents the power-transfer curve when the input power at the two turning points is similar, so the multivalued interval is very small.Fig. 7(a) presents the variation of the output power versus the coupling factor k for C p = 41.5 pF and different values of input power at the constant input frequency 15 MHz.The locus of turning points calculated in Section III-D has been superimposed.Fig. 7(b) presents the same analysis at the input frequency 14 MHz, with a smaller size of the locus of turning points, exhibiting a cusp at a lower k.As gathered from the above results, both the input frequency and input power have a significant influence on the shape of the nonlinear-resonance curves, strongly affected by the location of the turning points.Thus, the capability to directly obtain the locus of turning points will be valuable for an efficient prediction of the overall circuit behavior.This will be the object of Section III-D.

D. Mechanism for the Turning Point
To obtain the turning points, we will split (22) into real and imaginary parts and calculate the Jacobian matrix of the resulting two real equations with respect to |V 2 | and φ.This matrix should be singular at the turning points, and the singularity condition is formally identical to the one in (3) where the superscripts indicate real and imaginary parts.
As done in Section II, to get insight into the turning point mechanism, the charge Q(|V 2 |) will be decomposed into a linear term and a purely nonlinear contribution, that is, where the frequency dependent terms T 1 (ω), T 2 (ω), T 3 (ω) are And the following frequency-dependent components have been defined: Note that (25) has the same structure as (6).From the inspection of ( 26) and ( 27), the term T 1 (ω), associated with the linear contribution C|V 2 |, and the term T 2 (ω), affecting Q NL (d Q NL /d|V 2 |), are always positive.Under an equal sign of Q NL and d Q NL /d|V 2 |, the only possible negative contribution [capable to provide roots of (25)] is the one associated with ).This situation is analogous to the one obtained in the case of the single nonlinear resonator of Section II, where the term affecting ) is the only one that can lead the determinant in (7) to a zero value.For a more in-depth analysis of T 3 (ω), the common case of two equal inductors L p = L s = L will be assumed, so this term becomes where The frequency variation of T 3 for different k values is shown in Fig. 8.For low ω, the first and third terms of T 3 (ω) are negative and will dominate the whole expression.For high ω, the second positive term will dominate.If both Q NL and d Q NL /d|V 2 | are positive (negative) for all |V 2 |, one can expect to obtain turning points at lower (higher) frequencies than ω o .However, one must also note that, from certain k, there are three crossings through zero of T 3 (ω) and four sections with different signs, which can give rise to more than one turning point locus.For instance, in the case considered here having Q NL (|V 2 |) < 0 and ∂ Q NL (|V 2 |)/∂|V 2 | < 0, the turning point condition [singularity of (22)] requires T 3 (ω) > 0 which is obtained in two frequency intervals.As a result, there can be a coexistence of turning point loci that will be demonstrated here both theoretically and experimentally.
Using the same charge model And solving for |V 2 | 2 , one obtains the turning point voltage To obtain |V 2,T (k, ω)| 2 > 0, the elements T 3 and Q 3 must have opposite signs and the radicand must be positive or equal to zero (at the cusp), which requires a sufficiently large magnitude of T 3 in comparison with T 1 and T 2 .In this situation, for given k and ω, two turning points are obtained, except at the cusp, where the radicand is zero.In the case Q 3 < 0 considered here, the negative sign of the square root in (30) provides the turning points in the upper section of the nonlinear-resonance curves, which, when increasing the input frequency, give rise to a downward jump to the lower section.In turn, the positive sign provides the turning points in the lower section, which, when reducing the input frequency, give rise to an upward jump to the upper section.
Note that for each k (or ω) there can be one or more ω (or k) fulfilling the turning point condition.As in the case of the single resonator, the voltage amplitude at the turning points is scaled with √ |Q 3 |.Note that each |V 2,T | resulting from (30) corresponds to a different input power, which is directly obtained by replacing |V 2,T | in (23).
The turning point loci provided by (30) has been superimposed on the curve families in Fig. 6(a) and (b), providing the output power versus ω for two capacitor values, and in Fig. 7(a) and (b), providing the output power versus k for two ω values.They accurately pass through all the tuning points of the solution curves, corresponding to different P in values.
To obtain the locus of turning points in the plane defined by k, ω, one should perform a double sweep in these two parameters, calculate |V 2,T (k, ω)| from (30) and P in,T from (23).This will produce one or more surfaces of turning points P in,T (|V 2,T |, ω)in the space defined by k, ω, and P in .To summarize, the surfaces P in,T (|V 2,T |, ω) are obtained through the following three-step calculation: Note that though the parameters ω, k have been considered in (31), the surfaces can equally be calculated in terms of any two arbitrary parameters η 1 and η 2 .
When applying the new method to the circuit in Fig. 4, considering η 1 = k and η 1 = ω, and one obtains the two surfaces in Fig. 9(a).From these turning point surfaces one can obtain the family of turning point loci resulting from the variation of the input power.In fact, the locus of turning points for a particular input power, denoted as P in,o , is given by the contour plot of the surface(s) corresponding to P in,o .This can be seen in Fig. 9(b), which shows the family of turning point loci (obtained for different P in values) in the plane defined by ω and k resulting from the surfaces in Fig. 9(a).A relevant advantage of this method over continuation-based ones used in in-house software [14]- [17] is that it can predict disconnected loci sections.This is because it exhaustively provides all the loci points contained in the exploration intervals considered when calculating the surface.
As previously explained, the existence of two distinct turning point loci in the circuit of Fig. 4 loci can be understood from the form of variation of the term T 3 (k, ω), exhibiting four sign changes as shown in Fig. 8.The upper sections of the loci in Fig. 9(b), obtained for k > 0.5, correspond to a very small distance between the inductors.Thus, we will focus on the lower sections, with an expanded view in Fig. 9(c).The loci points are in total consistency with the curves shown in Figs. 6 and 7.This can be gathered by comparing, for instance, the frequency values at the turning points in Fig. 6(b) with the locus predictions for k = 0.15 and P in = 7.4, 10.96, and 13.52 dBm.A horizontal axis at k = 0.15 has been traced to facilitate the comparison.One can also compare the turning points in Fig. 7(a) and (b), with the locus predictions at f = 14 and 15 MHz) for the power values P in = 7.4, 10.96, and 13.52 dBm.Two vertical axes at the respective frequencies have been traced to facilitate the comparison.From the inspection of the loci family, one can understand the extremely small multivalued interval obtained for k = 0.15 and 14.99 MHz in Fig. 6(d).Note that the turning points in Figs. 6 and 7 had been validated with HB, which also constitutes a validation of the results in Fig. 9.
The loci provide at a glance an entire characterization of the multivalued regions of the nonlinear-resonance curves and thus, valuable global information on the geometry of these curves.Thus, they provide the evolution (versus any relevant IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES parameters) of the first turning point that gives rise to the upward jump from which the flattening is observed.From the inspection of Fig. 9(c), at a given input power, to obtain this turning point at a lower k (and, thus, a larger distance between the coils), one should reduce the input frequency.This agrees with the results of Fig. 7.As seen in Fig. 9(c), the turning points are obtained for lower k when using a higher input power.One can also gather that the multivalued region (comprised between the lower and upper point of the locus) decreases with the excitation frequency and increases with the input power.
As stated, at the cusp points the radicand of (30) is equal to zero, so they are given by the condition The voltage amplitudes are calculated from and scale with √ |Q 3 |.Condition (32) defines one or more curves in the plane k c , ω c , which are independent of the nonlinearity.They have been represented in Fig. 9 and accurately pass through the cusps of all the turning point loci.When approaching the cusp point the multivalued interval decreases, in consistency with Fig. 6(b) and (c), and as also gathered from the comparison of Fig. 7(a) and (b).As seen in Fig. 9(c) to obtain the cusp point at low k one needs a higher input power.On the other hand, as the input power decreases the locus becomes narrower and approaches the envelope of cusp points (cusp trend).

IV. NEW HB METHOD AND EXPERIMENTAL CHARACTERIZATION OF THE NONLINEAR-RESONATOR CIRCUIT
The new HB method relies on the same concepts used in the analytical studies of Sections II and III.Instead of using ordinary sweeps (as in a default simulation in commercial HB) or parameter-switching techniques (for instance, when optimizing an auxiliary generator [17]), the new methods are based on the calculation of the nonlinear function H (V, ω).This provides the response to a voltage amplitude V at the frequency ω at the analysis node.With HB, the function H (V, ω) will be calculated in a realistic manner, using detailed models of the circuit elements (including parasitic effects) and as many harmonic components as desired.The results of the new HB method will be validated with default HB simulations and through an experimental characterization of the linear resonator coupled to a nonlinear one.

A. New HB Method
As stated, the new HB method is based on the calculation of H (V, ω), instead of directly addressing the circuit with its input source included.Once the function H (V, ω) is available, one will obtain the solution curves by means of the complex equation E g e j φ = H (V, ω), in a manner like what was done in Section III.To obtain the function H (V, ω) in commercial HB one must carry out a single (simultaneous) simulation of Fig. 10.Sketch of the analysis procedure to obtain the solution curves and turning point loci of the linear resonator coupled to a nonlinear one.Note that the two circuits are analyzed in a single simulation.(a) Circuit with E g = 0 is used to calculate the total nonlinear admittance function Y (V, ω) at Node 2. (b) Circuit is used to obtain the linear function F(ω) that relates the Norton equivalent current of the input network with the input source E g .two circuits.The sketch of this simulation, presenting the two circuits that are analyzed at once, is shown in Fig. 10.In the first circuit, we set the input source to zero and make use, instead, of an auxiliary generator.This is a voltage source with the amplitude V at the frequency ω, in series with an ideal bandpass filter at ω [17].The auxiliary generator is connected in parallel between the upper node of the antiseries diodes (Node 2) and ground, so the amplitude V agrees with |V 2 |.The circuit with the auxiliary generator is used to obtain the total nonlinear admittance function Y (V, ω) at Node 2, so V = |V 2 |.The function Y (V, ω) is calculated through a double sweep in ω and V , with the pure HB system acting as an inner tier.The second circuit, also shown in Fig. 10, is used to obtain the Norton equivalent of the linear input network at Node 2. A voltage source with the dummy value E g1 = 1 V is introduced in this linear circuit, and a function F(ω) is calculated through the relationship The two circuits in Fig. 10 are simulated simultaneously to obtain Y (V, ω), from the first circuit, and F(ω), from the second circuit.Once the two functions Y (V, ω) and F(ω) are available one can calculate the desired nonlinear function H (V, ω).To obtain this function, one must consider that the circuit solutions should fulfill the combined system Y (V, ω)V e − j φ = I e j ϕ I e j ϕ = F(ω)E g (35) where I e j ϕ is the complex Norton-equivalent current at the analysis node and −φ is the phase of the node voltage.The two relationships in (35) can be compacted in the following single complex equation: Thus, by performing a double sweep in ω and V , and carrying out the simultaneous HB simulation of the two circuits in Fig. 10, one can obtain H (V, ω) as Note that in this calculation the number NH of harmonic terms can be as high as desired, since Y (V, ω) is, in fact, an outer-tier function: for each pair of excitation values ω and V , the pure HB system acts as an inner tier.Once the function H (V, ω) is available, one can obtain all the solution curves V versus ω at arbitrary E g values from the corresponding contour levels of the surface |H (V, ω)|.On the other hand, the locus of turning points is given the contour level of zero value of ∂|H (V, ω)| 2 /∂ V , easily obtained in the output window of the commercial software using a differentiation function.If, instead of ω, a different parameter η is used, one should perform a double in sweep η and V to calculate the function H (V, η) and then obtain the solution curves in terms of V versus η at different E g values from the corresponding contour levels of the surface |H (V, η)|.
To calculate the turning point locus in the plane defined by any two arbitrary parameters η 1 , η 2 , we will apply the same three-step procedure of (31).We will perform a sweep in η 1 and, at each point of the sweep, we will obtain the turning point locus in terms of V and η 2 , given by the condition: ∂|H (V, η 2 )| 2 /∂ V = 0. To translate this turning point locus to the plane defined by η 2 , P in , one should calculate the input power at each locus point V T , η 2 , by applying The collection of turning point loci resulting from the sweep in η 1 provides a turning point surface in the space η 1 , η 2 , P in .Then the turning point locus at any desired P in value is directly obtained from the corresponding contour level of the turning point surface.

B. Simulations and Measurement Results
The design of the coupled inductors has been carried out using Wheeler's equation [23], [24].The prototypes have been manufactured on FR4 (h = 62 mil) substrate, while the supports for the coils have been machined through computer numerical control (CNC) on extruded acrylic clear plastic sheets (4 mm thickness).Several inductor pairs were manufactured to ensure the actual inductance values were closer to the ones used in the circuit simulation.The final implementation parameters are given in Table I.
The measurement set up is shown in Fig. 11.A WW2572 waveform generator (Tabor Electronics, Nesher, Israel) has been used as the input source.The voltage at the output node has been measured using a DSO6034A  11.Experimental setup.The primary and secondary resonators have been manufactured using FR4 (h = 62 mil) substrate.Coils have been mounted on supports machined from acrylic clear plastic sheets using AWG 18 copper wire.A WW2572 waveform generator (Tabor Electronics) has been connected to the linear resonator and the voltage waveforms have been measured using a DSO6034A oscilloscope (Agilent Technologies.)oscilloscope (Agilent Technologies, Santa Clara, CA, USA).For a more accurate acquisition, the output voltage waveforms are averaged, and, for the channel measuring the voltage at the load resistor, the BW limit option of the digital storage oscilloscope (DSO) has been enabled.This allows a more accurate estimation of the voltage at the fundamental frequency while the higher harmonic terms are filtered.The automated acquisition measurement is primarily based on an increasing input-voltage sweep for each test frequency and fixed separation of the coupled coils.When a discontinuity is detected due to a jump caused by the lower-section turning point, a decreasing sweep in the input voltage is also carried out to complete the upper section of the solution curves.

C. Validation With Experimental Results
In the following analyses, the varactor diodes SMV1470-004LF are described with their full manufacturer models, including all the parasitic elements.When considering these full diode models, the circuit element values have been slightly varied for an optimized response.The most relevant modifications are the change of the output resistor to R L = 2.5 k and the capacitor C p to 39 pF.All the simulations through the new method are carried out using NH = 7 harmonic terms.Fig. 12(a) presents the power-transfer curves obtained for k = 0.15 at the input frequencies 14, 15, 16, and 17 MHz.The results are compared with those provided by the analytical formulation of Section III.The charge polynomial model has been slightly modified for a better match in the presence of full diode models and is now given by C = 55 pF and Q 3 = −5.8C/V 3 .As can be seen, there is quite a good agreement.Fig. 12(b), showing the same simulation with default HB in commercial software, evidences the serious convergence difficulties of this simulation.There is a large input-voltage interval with no convergence (and, therefore, no solutions), which prevents the accurate prediction of the power-transfer curve Fig. 12(c)-(f) compared with the power-transfer curves obtained through the numerical method with experimental measurements at the respective input frequencies 14, 15, 16, and 17 MHz.For these measurements, which are based on the setup of Fig. 11, the coupling factor has been experimentally characterized versus the distance between the coils.From this characterization, the value considered in simulation, k = 0.15, approximately corresponds to the distance d = 7.1 cm between the coils.In Fig. 12(c)-(f), the accuracy of the new method is also validated with default HB simulations (unable to pass though the turning points).
Fig. 13(a) presents the variation of the voltage amplitude |V 2 | versus the input frequency for different input powers at k = 0.15.Experimental measurements are superimposed.As can be seen, there are two distinct turning point loci, giving rise to two different folding effects, at lower frequencies (though higher than ω o ), for the smaller P in values, and at Measurements are superimposed.(c) Exhaustive experimental characterization of the circuit response versus the input frequency when considering the input power interval P in = −2.26dBm to P in = 7.5 dBm.higher frequencies, for the higher P in values.This is explained from the resonance analysis in Section III-B.As the input power increases, the first resonance folds rightwards, then it is the second resonance the one that folds rightwards.
For better insight into the mechanisms giving rise to the two loci, the same analysis has been carried out with the analytical model, which provides the results in Fig. 13(b), with experimental measurements superimposed.As can be seen, two loci are also obtained with the analytical formulation, in which the back-to-back diodes are described with the same cubic-polynomial model considered in the single resonator analysis of Section II.Thus, the existence of the two loci is  due to the frequency dependence of linear circuitry, which is more complex than the one in a single nonlinear resonator.This was anticipated in Section III from the analysis of the frequency variation of the term T 3 (ω) of the singularity condition (25), shown in Fig. 8. Fig. 13(c) presents an exhaustive experimental characterization of the circuit response versus the input frequency when considering the input power interval P in = −2.26dBm to P in = 7.5 dBm.This characterization has been carried out through a fine frequency sweep.The envelopes of the jump points (experimental turning point loci) are superimposed.
Several factors may cause the discrepancies between the simulation and experimental results.Tests were carried out considering variations of the package parasitics within reasonable tolerance limits and we observed that slight variations gave rise to nonnegligible discrepancies between the resulting solution curves.On the other hand, the spiral inductors were manufactured following Wheeler's equations [23] for planar coils.As shown in [24], these equations have an estimation error up to 20%, depending on different factors (e.g., number of turns, spacing, etc.).The measured self-resonance frequency of the inductors is about 30 MHz, and parasitic effects were present in the operation frequency range (11)(12)(13)(14)(15)(16)(17).
Fig. 14 presents the efficiency curve versus the coupling factor k (lower axis) at the input frequency f in = 15 MHz and power P in = 5.1 dBm, where the results are compared with those obtained with the analytical formulation.Experimental measurements versus the distance d between the coils for the same input frequency and several P in values are also shown.To avoid the impact of any inaccuracies in the estimation of the relationship between the coupling factor and the distance, an additional upper axis indicating this distance has been introduced in Fig. 14.Note that the objective of this work is not to obtain a demonstrator of a very performant power-transfer system, but to derive new simulation methodologies that should facilitate the design and optimization of these systems.Thus, no effort has been devoted to optimizing the inductor coupling and or modifying the nonlinear capacitance to maximize the flattening of the response.
Finally, Fig. 15 presents the turning point loci obtained for P in values comprised between 1.5 and 6 dBm through a previous calculation of a turning point surface.They accurately predict the two different folding effects at lower input frequency (for lower input power) and at higher input frequency (for higher input power).A horizontal axis at k = 0.15 has been introduced to facilitate the comparison of the loci predictions with the turning points in the curves of Fig. 13(a) traced versus frequency for different P in values.A vertical axis at 15 MHz has also been introduced to facilitate the comparison of the loci predictions with the turning points in the numerical curve of Fig. 14 traced versus k.

V. CONCLUSION
A detailed frequency-domain analysis of circuits based on a nonlinear resonator has been presented, illustrated through its application to two configurations recently proposed in the literature for energy harvesting and wireless power transfer.In the two cases and considering a describing-function representation of the nonlinear charge, explicit expressions of the turning points that give rise to the folding of the resonance curves have been derived.These expressions provide insight into the turning point mechanism and its dependence on the circuit elements and parameters.The mathematical condition fulfilled by the cusp points, from which the solution curves are no longer multivalued, has also been derived.For a realistic prediction of the behavior, a new single-step method to obtain the multivalued solution curves in commercial HB has been provided, together with a new procedure to determine the locus of turning points in terms of any two analysis parameters.This is done through the calculation of a turning point surface in the space defined by the two parameters and the input power.The methods have been applied to a practical wireless powertransfer system at 15 MHz, obtaining a very good agreement with the experimental results.
Fig.3.Analysis of the nonlinear resonator in Fig.2.(a) Nonlinear resonance curves for different values of the input-current amplitude, in terms of V versus ω, for G = 1/500 = 2 mS.Default HB simulations are superimposed for I g = 9 mA.HB results obtained with the numerical method in Section IV are also presented for comparison.The locus of turning points obtained through (9) is superimposed and the cusp resulting from (11) and (12) is indicated.(b) Same for G = 1/10 3 = 1 mS.Default HB simulations are superimposed for I g = 4.5 mA.(c) Evolution of the turning locus when varying the conductance G.The dashed line provides the evolution of the cusp.
Authorized licensed use limited to: BIBLIOTECA DE LA UNIVERSIDAD DE CANTABRIA.Downloaded on July 29,2021 at 09:02:58 UTC from IEEE Xplore.Restrictions apply.

Fig. 5 .
Fig. 5. Linear operation of the two coupled resonators in Fig. 4. (a) Evolution of the circuit poles in the complex plane for L p = L s = 3.2 μH and C p = C = 55 pF.(b) Evolution of |H L ( jω)| versus k for the same values.(c) Variation of the imaginary part of the poles versus k for the same values and for L p = L s = L = 1.6 μH and C p = C = 83 pF.(d) Evolution of |H L ( jω)| versus k for L p = L s = L = 3.2 μH, C p = 41.5 pF, and C = 55 pF.

Fig. 6 (
a) and (b) presents the family of output-power curves of the circuit in Fig. 4, obtained versus ω for different values of input power P in = E 2 g /(8R s ), at the constant coupling coefficient k = 0.15 for L p = L s = 3.2 μH and C p = C = 55 pF, in Fig. 6(a), and for L p = L s = 3.2 μH, C p = 41.5 pF and C = 55 pF, in Fig. 6(b), which are the two sets of circuit-element values considered in Section II-B.The results are compared with HB simulations, through the numerical method in Section IV, based on the calculation of |H (|V 2 |, ω)| in commercial HB, plus the tracing of contour levels.The locus of turning points calculated in Section III-D is also superimposed.Note the significant BW enhancement and flattening of the resonance curves in comparison with the linear case of Fig. 5(c) for the same k value.This is due to the broadening of the lower-frequency resonance and the folding of the upper-frequency one.The solution curve at a constant input frequency ω is obtained by simply sweeping |V 2 | and calculating E g from the absolute value |H (|V 2 |, ω)| [see(22)].This has been done

Fig. 6 .
Fig. 6.Linear resonator coupled to a nonlinear one with a simplified diode model.The coupling factor is k = 0.15.(a) C p = 55 pF.Family of nonlinear resonance curves obtained for input power comprised between 4.39 and 13.5 dBm.The turning point locus is superimposed.The analytical results are compared with HB simulations using the method in Section IV.(b) Same for C p = 41.5 pF.(c) Power-transfer curves for C p = 41.5 pF considering three different input frequencies.The default HB is unable to pass through the turning points.Simulations using the numerical method in Section IV are superimposed.(d) Power-transfer curve for C p = 41.5 pF and f = 14.99 MHz, with a very small multivalued region.

Fig. 7 .
Fig. 7. Linear resonator coupled to a nonlinear one with a simplified diode model.Variation of the output power versus the coupling factor k for C p = 41.5 pF and different values of the input power.The turning point loci obtained with the method described in Section III-D are superimposed.(a) At the constant input frequency 15 MHz.(b) At 14 MHz.

Fig. 9 .
Fig. 9. Calculation of the turning point locus in terms of ω and k for different P in values.(a) Newly defined turning point surfaces in terms of the analysis parameters ω, k, and P in .The envelope of cusp points defined by (32) is shown in dash-dotted line.(b) Family of turning point loci in the plane defined by ω and k for different input powers.The envelope of cusp points is superimposed.(c) Expanded view of the turning point loci obtained for the lower k values.

Fig. 12 .
Fig. 12. Linear resonator coupled to a nonlinear one when describing the varactor diodes with realistic models.Power-transfer curve for different values of input frequency at k = 0.15 (d = 7.1 cm).(a) Comparison with the results obtained with the analytical formulation using the charge model Q(V ) = C V + Q 3 V 3 .(b) Default HB simulation in commercial software, evidencing serious convergence difficulties.(c) Simulation with the new method and default HB with measurements superimposed at the input frequency 14 MHz.(d) Same for 15 MHz.(e) Same for 16 MHz.(f) Same for 17 MHz.

Fig. 13 .
Fig. 13.resonator coupled to a one when describing the varactor diodes with realistic models.Variation of the voltage amplitude |V 2 | versus the input frequency at k = 0.15 (d = 7.1 cm), when varying the input power.(a) New HB-based method with accurate diode models.There are two turning point loci.Measurements are superimposed.(b) Analytical formulation, when describing the diodes with the model Q(V ) = C V + Q 3 V 3 .Measurements are superimposed.(c) Exhaustive experimental characterization of the circuit response versus the input frequency when considering the input power interval P in = −2.26dBm to P in = 7.5 dBm.
Authorized licensed use limited to: BIBLIOTECA DE LA UNIVERSIDAD DE CANTABRIA.Downloaded on July 29,2021 at 09:02:58 UTC from IEEE Xplore.Restrictions apply.

Fig. 14 .
Fig.14.Linear resonator coupled to a nonlinear one when describing the varactor diodes with realistic models.Variation of the power-transfer efficiency E ff versus the coupling factor k (upper axis) for 15 MHz and P in = 5.1 dBm.The results are compared with those obtained with the analytical formulation.Experimental curves versus the distance d between the coils (lower axis) are also presented.

Fig. 15 .
Fig.15.Linear resonator coupled to a nonlinear one when describing the varactor diodes with realistic models.The analysis parameters are η 1 = k, η 2 = ω.Family of turning point loci obtained for input powers comprised between 1.5 and 6 dBm.It accurately predicts the two different folding effects.

TABLE I DESIGN
PARAMETERS FOR COUPLED INDUCTORS Fig.