Nonlinear Analysis of Oscillator Mutual Injection Locking Through Inductor Coupling

This work presents an in-depth investigation of the nonlinear behavior of two mutually injection-locked oscillators through inductor coupling. An analytical formulation, solved through an innovative procedure, facilitates the understanding of the qualitative transformations in the system solutions when increasing the coupling factor. The analysis demonstrates that, in a manner similar to the unilaterally injection-locked oscillators, families of disconnected/connected curves are obtained when increasing this factor, although the patterns, associated with distinct operation modes, are more complex. Then, an accurate numerical method to predict the behavior of coupled transistor-based oscillators is presented, based on nonlinear admittance models of the individual oscillators. Mathematical conditions are derived to solve the coupled system through a two-level contour-intersection technique. In this way, all the solutions coexisting for a given set of element and parameter values are calculated simultaneously, in an exhaustive manner. The cases of two coupled oscillators at the fundamental frequency and at 1:3 frequency ratio are considered. Possible applications include the oscillator phase-noise reduction and the implementation of sensors using the phase shift between the two oscillator elements.


I. INTRODUCTION
R ECENTLY, successful implementations of mutually injection-locked oscillators through inductor coupling have been demonstrated [1]- [14]. In a two-oscillator system, the tuning range versus a tuning parameter is increased due to the two resonance frequencies exhibited by the coupled resonators [7]- [9]. Other applications include sensors [4], [10]- [13] and near-field wireless data systems in the vein of [14] at 13.35 and 30 MHz. In addition, the oscillator coupling enables a reduction of phase noise, which follows the rule S N = S − 10log 10 N [15] where S is the spectral density of a single oscillator [16] and N is the number of oscillator elements.
Despite the interest in inductively coupled oscillators, there is a lack of insightful analysis tools for these configurations, in which several oscillation modes may coexist. Murphy and Darabi [3] propose a useful eigenvector/ eigenvalue analysis [17]- [19] of the full oscillator system, which takes advantage of the circulant structure of the coupling matrix. The coexistent oscillation modes are identified; however, the analysis described is not valid in the presence of asymmetries. In turn, Suárez et al. [20] present an investigation of two mutually injection-locked oscillators that is not restricted to identical oscillators. However, the analysis relies on a linearization of the total admittance functions of the individual oscillators about their stand-alone free-running solutions [21]- [26], which limits its accuracy to weak coupling conditions [27], [28].
As shown in this work, nonlinear effects can be observed from rather small values of the inductive coupling factor k. To account for these effects, we will extend [20] by addressing the coupled system under significant deviations of the oscillator elements with respect to their free-running solutions. Initially, an analytical formulation is derived, which provides insight into the qualitative transformations undergone by the system when increasing the mutual coupling. Each oscillator is affected by the other, so the behavior is more complex than in oscillators injection locked by an independent source [29], [30]. As will be demonstrated, families of disconnected curves that merge from a certain k value are obtained versus a tuning parameter.
To analyze a realistic system of transistor-based coupled oscillators, one must keep in mind the coexistence of solutions revealed by the analytical investigation. In fact, unveiling the behavior pattern of inductively coupled oscillators is the main contribution of this work. The use of time-domain integration would require a global exploration of all the possible initial values (which is virtually impossible) since each stable solution has its own basin of attraction [31], [32]. The problem would be similar in the envelope domain [33]- [35]. On the other hand, when using harmonic balance (HB), one must be able to lead each oscillator to an oscillatory state, which can be done connecting an auxiliary generator (AGs) to each oscillator element [29]- [36]. The amplitude of each AG, together with their phase shift and the oscillation frequency, must be optimized to fulfill (simultaneously) the oscillation condition 0018-9480 © 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
at the two oscillator elements. However, this procedure often fails due to the need to fulfill four goals, corresponding to the zero value of the real and imaginary part of the admittance function in each oscillator [18], [29], using only two observation/analysis nodes. Here, an alternative numerical method will be presented. The oscillators are described with numerical nonlinear admittance functions, extracted from HB, while the passive coupling network is described with a linear admittance matrix. Advantage is taken of the fact that the oscillator admittance functions do not depend on the excitation phase. Instead, the phase shift between the oscillator elements is considered in the formulation of the entire system, which is solved separately from HB, using custom software. The methodology constitutes the first generalization of the semianalytical formulation of coupled-oscillator systems [21], [22] to oscillators described with nonlinear admittance functions. Furthermore, judicious mathematical conditions allow extending the contour-intersection technique in [37] and [38] to the coupled system. This procedure (which is not based on continuation [39]- [41]) allows the simultaneous and exhaustive detection of all the oscillation modes coexisting for a given parameter value. This article is organized as follows. Section II summarizes the linearized analysis of two inductively coupled injectionlocked oscillators at the ratio 1:1 presented in [20]. Section III describes the nonlinear analytical investigation of this system. Section IV presents the numerical analysis method, applicable to realistic transistor-based oscillators. Section V describes the extension of this method to the oscillator coupling at the ratio 1:N.

II. TWO MUTUALLY INJECTION-LOCKED OSCILLATORS
AT 1:1 WITH LINEAR OSCILLATOR MODELS A. System Description Let a general system of two inductively coupled oscillators at the frequency ratio 1:1 be considered. The system is shown in Fig. 1(a). The passive-linear coupling network (in dashed line) will be described with its admittance matrix where L 1 and L 2 are the inductors of the first and second oscillator and M = k·(L 1 L 2 ) 1/2 is the coupling inductance. Note that, in the general case, the admittance matrix y i j may include parasitics. The oscillators are described by their individual admittance functions Y T 1 (V, ω) and Y T 2 (V, ω), where V and ω are the amplitude and frequency of the excitation voltage. Assuming mutually injection-locked operation at the frequency ω, the coupled system can be described (at the fundamental frequency) with the complex system where V 1 and V 2 are the oscillation amplitudes and φ is the phase shift between the two voltage signals. Note that the admittance − j /(L i ω) associated with the coupled inductor in each oscillator has to be subtracted from the functions Y T 1 and Y T 2 since this inductive effect is included in (1).
Several previous works [19]- [26] have successfully developed semianalytical formulations of complex systems based on admittance-type oscillator models extracted from HB simulations. However, in all these works, the admittance-type models are linearized about their corresponding stand-alone free-running solutions, which is valid under the assumption of a small deviation from these solutions once each oscillator is introduced into the system. In the case of (2), the two admittance functions Y T 1 (V 1 , ω) and Y T 2 (V 2 , ω) are linearized [20] about their respective free-running solutions, which, as will be shown, limits the validity of the analysis to small k values. Here, we will address for the first time the nonlinear operation of the inductively coupled system, considering the two bivariable nonlinear functions Y T 1 (V 1 , ω) and Y T 2 (V 2 , ω), instead of their linearized approximations. To facilitate the comparison between the linear and nonlinear cases, Section II-B summarizes the results obtained in [14], with linearized oscillator models.

B. Formulation Based on Linearized Models
The system based on linearized oscillator models (valid for small k) will allow obtaining the variation of the coupledsystem solution versus a parameter η (a capacitor value, for instance). The frequency and amplitude of the free-running solution (in stand-alone operation) of each of the two distinct oscillators are given by ω oi and V oi , where i = 1, 2. In practice, these free-running solutions can be calculated with the aid of an AG [ Fig. 1(b)], optimizing its amplitude V AG and frequency ω AG to obtain a zero value of the AG current-to-voltage: Y AGi , where i = 1, 2 [18], [29]. In fact, the AG admittance function agrees with the defined oscillator-admittance function Y T i (V i , ω). In the linearized analysis, this function is described with its first-order Taylor series expansion about the freerunning solution fulfilling Y T i (V oi , ω oi ) = 0. The amplitude derivative ∂Y T i /∂ V i is calculated by setting the AG frequency to ω AG = ω oi and considering a small increment in the voltage amplitude V AG about V oi [17]. Note that a full HB analysis, with as many harmonic components as desired, is carried out after the application of the increment in V AG . Likewise, the frequency derivative ∂Y T i /∂ω is calculated by setting the AG amplitude to V AG = V oi and considering a small increment in the AG frequency about ω oi [23], [24]. Finally, the derivative ∂Y T i /∂η will be obtained through the same procedure: considering a small increment in η about the freerunning point V oi , ω oi .
Synchronized operation of the inductively coupled oscillator system at the frequency ω will be assumed. Under a small variation of the parameter η, it will be approximately described as follows: where for simplicity the two oscillators are assumed to have the same free-running frequency ω o1 = ω o2 . In addition, the following quantities (evaluated at ω o ) have been defined: ∂Y To obtain the solution, (3) is split into real and imaginary parts, which provides a system of four real equations in four unknowns: V 1 , V 2 , φ, and ω. In practice, the system is solved by sweeping φ, and solving for V 1 , V 2 , ω, and η.
The above analysis has been applied to two oscillators of the Van der Pol type [shown in Fig. 1(c)], with the element values a = −0.03 A/V, b = 0.01 A/V 3 , L 1 = L 2 = 33 nH, C 2 = 76.7 pF, R 1 = R 2 = 100 . The two oscillators are assumed identical, except for C 1 = C 2 + C 1 , acting like the analysis parameter η = C 1 . System (3) only enables a valid prediction of the steady-state solutions for very small k values. In the analysis of Fig. 2(a), the coupling factor is set to k = 0.01. In Fig. 2(a), the phase φ has been represented versus C 1 . The results of the analytical formulation are compared with an HB simulation with 15 harmonics (connecting an AG to each oscillator). The excellent agreement is because the derivatives of Y T 1 and Y T 2 are calculated with the HB system as an inner tier [26], considering the same number of harmonic terms. As will be shown in Section III, under a coupling factor as small as k = 0.1, the coupled system The stability and phase-noise analyses are based on the introduction of a small perturbation of complex frequency in system (3), which involves a subsequent linearization with respect to each solution obtained versus C 1 , as shown in [19].

III. NONLINEAR ANALYTICAL FORMULATION OF TWO MUTUALLY INJECTION-LOCKED OSCILLATORS
The inductive coupling factor can take any value 0 < k < 1, so in most cases, it will not be possible to linearize the admittance functions Y T 1 (V, ω) and Y T 2 (V, ω) of the two individual oscillators about their respective free-running solutions. In the following derivations, this linearization is avoided, and the oscillators are described with nonlinear admittance models.

A. Formulation
The nonlinear behavior of the two inductively coupled oscillators can be understood by particularizing (2) to the system of two Van der Pol oscillators in Fig. 1(c). The cubic nonlinearity in each oscillator is modeled with its corresponding describing function limited to the fundamental frequency, which provides the following matrix system: where L 1,2 = L, R 1,2 = R, G T = a + 1/R and, for compactness, the following definitions have been introduced: The complex system (5) contains four unknowns V 1 , V 2 , ω, and φ. The initial objective is to eliminate the phase shift φ and derive a system of three real equations in the remaining three unknowns V 1 , V 2 , and ω. It is taken into account that (5) is homogeneous, so its associated matrix must be singular, which provides where the following additional definition has been used: From the first row of (5), one can derive the complex equation Now, splitting the complex equation (7) into real and imaginary parts and obtaining the squared magnitude of (8), one obtains the following system of three real equations in the three real unknowns V 1 , V 2 , and ω: In the above system, one can make the amplitudes V 1 and V 2 disappear and, thus, obtain a real equation in the frequency ω The above equation cannot be solved explicitly but provides the two error functions (H + and H − ), respectively, corresponding to the plus and minus signs before the square root. The zeroes of (10) provide the frequencies of the potential solutions. Note that additional conditions on V 1 and V 2 must be fulfilled. The error functions H ± obtained for C 1 = 77 pF and several values of the coupling factor k are shown in Fig. 3. For each k, the frequency values of the potential solutions correspond to the crossings through zero.
Once the potential-solution frequencies are known, one can apply a straightforward procedure to obtain the solution amplitudes. Using (9)(a) and (9)(b), one derives the following directly solvable equations: On the other hand, actual solutions must also fulfill (9)(c) Once the amplitudes V 1 and V 2 are known, the phase shift between the two elements is obtained from (8).
In case the two oscillators are identical, C = 0, one obtains the two in-phase (0 • ) and out-of-phase (180 • ) modes with the identical amplitudes V 1 = V 2 resulting from the eigenvalue analysis in [3] and [17]. However, as will be shown here, and for relatively low values of the coupling factor k, there will also be two other modes with different amplitudes (V 1 = V 2 ) in the two oscillators. The frequency of the nonsymmetric modes is determined by the condition which provides a particular null of the error function (10). From (13), the frequency of the nonsymmetric solutions is The phase shift values can be obtained from the first equation in (5) where (13) has been taken into account. Because the amplitudes V 1 and V 2 are positive and different, there can only be two possible values of phase shift, given by φ = ±π/2. Thus, the amplitudes of the remaining solutions fulfill where the upper (lower) signs in (16)(a) and (16)(b) are associated and correspond to the same solution. Note that the relationships (16) can also be derived from (12). These nonsymmetrical solutions under identical oscillators are obtained here for the first time to our knowledge.

B. Behavior When Increasing the Coupling Factor k
Using the above procedure, we have obtained the solution curves of the two inductively coupled Van der Pol oscillators versus the capacitance C 1 (in the first oscillator) when increasing k. The potential-solution frequencies are calculated performing a double sweep in C 1 and ω, and obtaining the zero-value contours of H ± . Then, the solution curves are achieved using the additional equations (11) and (12). When considering different values of the coupling factor k = 0.1, 0.12, and 0.2, one obtains the evolution shown in Fig. 4 (k = 0.1), Fig. 5 (k = 0.15), and Fig. 6 (k = 0.2). In each case, the oscillation frequency, amplitudes, and phase shift are represented in (a), (b), and (c), respectively, versus C 1 . For validation, the results are compared with circuit-level HB simulations using two AGs. However, one must note the following: 1) the HB convergence is facilitated by the simple topology of the cubic-nonlinearity oscillators and 2) the HB solution curves could only be obtained by first providing the values resulting from the analytical formulation to the AGs used in this HB simulation. The AG optimization confirmed the validity of these analytical solutions. Then, the solution curves could be completed in HB due to the inherent continuation procedure of the HB sweep.
For the lower k values (Fig. 4), one obtains three disconnected solution curves. One is a closed eight-shaped curve, whereas the other two curves are open. Note that in the  Fig. 4(a)], the higher amplitude is V 1 (in the order of 1.7 V) and the lower one is V 2 (below 0.5 V). In the other open solution, the situation is opposite.
As gathered from Fig. 4(b), for each C 1 , there are three solution frequencies. This is can be seen at C = 0, where two of these frequencies agree with those of the two in-phase and out-of-phase modes ω = 1/(LC(1 ± k)) 1/2 (for which the amplitudes of the two oscillators are equal V 1 = V 2 ) [20], and the third frequency is ω = 1/(LC(1 − k 2 )) 1/2 (for which V 1 = V 2 ). Regarding the phase shift, and in agreement with (13)-(16), for two identical oscillators ( C = 0), besides the two modes with phase shifts 0 • and 180 • , there are two other modes with the phase shifts ±90 • , the same frequency ω = 1/(LC(1 − k 2 )) 1/2 and distinct amplitudes.
When increasing k, the closed eight-shaped curve splits into two sections and each section merges with one of the open curves [ Fig. 6], to give rise to two distinct open solution curves. In these open curves, each oscillator is dominant at each side of the middle value C = 0. The described curve merging can be compared with the structural behavior of a single oscillator injection locked by an independent source [29], [30], [42], [43]. In that case, when increasing the amplitude of the injection source, a single closed curve and a single open curve merge into a unique curve, which for the lower values of the input amplitude exhibits turning points or folding. Here, when increasing k, we obtain a double family of curves, one corresponding to each operation mode. In one of the two modes V 1 is higher (lower) than V 2 for C < 0( C > 0). The opposite is true of the other modes. The curves exhibit turning points for intermediate k values (Fig. 5), at which the HB simulations undergo a discontinuous jump. The curves become less intricate when increasing k, as gathered from Figs. 5 and 6.
The variation of the phase shift versus the analysis parameter C 1 is shown in Figs. 4(c), 5(c), and 6(c). As stated, for the lower k values and C = 0, there are always two modes with the respective phase shifts 0 • and 180 • , as well as two additional modes with ±90 • . In the closed solution curves versus C 1 , the phase variations are strong and cover the whole range −180 • to 180 • . When increasing k, the closed curves split and merge, and the modes with ±90 • at C = 0 disappear [ Fig. 5(c)]. In the open curves, the phase excursion decreases with k. Thus, for sensor applications, one should use a relatively small k.
To understand the stability properties of the coexisting modes, one should note that the eight-shaped solution curves obtained for the lower k values turning points (T 1 , T 2 , T 3 , T 4 ), as seen in Fig. 4. At each turning point, a real pole crosses through zero. Performing a stability analysis, one obtains that the section of the eight-shaped curve comprised between T 1 and T 4 is stable, and so is the section between T 2 and T 3 . Thus, there are two coexistent stable modes. At the turning points, the oscillators become unlocked and there is a transition to a doubly autonomous quasi-periodic regime. On the other hand, the open solution curves are always unstable.
After each eight-shaped curve splits and merges with one of the open solution curves, Hopf bifurcations will take place (from certain k) in the two single and distinct open curves. This is illustrated in Fig. 7(a), which presents the stability analysis of the two independent oscillation modes obtained for k = 0.2 using pole-zero identification [44], [45]. The real part of the dominant poles has been represented versus C 1 . The two modes are stable in the central interval about C = 0. This stable interval is bounded by secondary Hopf bifurcations [ Fig. 7(a)] at which a pair of complex-conjugate poles crosses the imaginary axis to the right-hand side of the complex plane [29], [42], [46]. When the Hopf bifurcation occurs, the system evolves into a self-oscillating mixer regime. The detected Hopf bifurcations have been superimposed in Fig. 6. The length of the stable intervals is different for the two modes. When the two modes are simultaneously stable, one or another will be experimentally observed depending on the initial conditions. This has been validated with the time-domain simulations in Fig. 7(b) and (c). Note that the coexistence of stable periodic modes with different phase shift is possible due to presence of unstable dc and quasi-periodic solutions [47], [48], acting as separators of the basins of attraction of the two coexistent stable modes.
To summarize, under small k values, there is a closed solution curve through which the two oscillators are in an oscillatory state and mutually injection locked, plus two other disconnected curves in which system behaves as if only the first (second) oscillator is in an oscillatory state and the second (first) oscillator is responding to this oscillation signal. From certain k, two distinct parts of the previously closed curve merge separately with the disconnected curves and give rise to two open curves at different frequencies, each corresponding to a different oscillation mode. After this merging, the two open curves exhibit strong folding due to the presence of turning points, or infinite-slope points. The higher phase sensitivity is obtained for the smaller k values, through the closed solution curves.

IV. NONLINEAR NUMERICAL ANALYSIS OF TWO MUTUALLY INJECTION-LOCKED OSCILLATORS
In this section, a general numerical analysis of two inductively coupled transistor-based oscillators will be presented. As explained below, for this numerical analysis, the coupled system is best formulated in the following manner: where the nonlinear-admittance functions Y A1 (V 1 , ω) and Y A2 (V 2 , ω) may not strictly correspond to the stand-alone oscillator circuits; they may represent just part of these oscillators. In that case, the passive admittance matrix [y p (ω, k, η)] accounting, in principle, for the coupling network (and depending on the coupling factor k), would also include part of the linear networks of these oscillators. For convenience, the analysis parameter η (a tuning capacitor, for instance) might be also included in [y p (ω, k, η)]. The subtraction of the linear elements to be included in the coupling network is performed after the calculation of those functions (with the aid of AGs).
The procedure will be illustrated through its application to the system of two differential bipolar-based oscillators shown in Fig. 8(a) and (b). When uncoupled, they oscillate at the frequency f o = 25 MHz (for C 1 = C 2 = 22 pF). To obtain each function Y Ai (V, ω), where i = 1, 2, each oscillator is simulated with an AG connected between the terminals at which this admittance function is defined, as shown in Fig. 8(c). Then, a double sweep in ω and V is carried out. The flexibility in the choice of the analysis terminals, as well as the order of the nested sweeps (with the one in V being the internal one), should avoid potential convergence problems of HB. In principle, system (17) could be solved through a Newton iteration, after a suitable analytical modeling of the nonlinear active admittance functions, though this is beyond the scope of this initial article. On the other hand, the advantage of the procedure presented in the following is that it enables an exhaustive search of solution curves and operation modes, which would not be possible through an ordinary error-minimization/continuation method [39]- [41].
To solve (17), advantage is taken of the fact that each of the admittance functions Y A1 (V 1 , ω) and Y A2 (V 2 , ω) depends only on the amplitude and frequency of its own excitation. The singular determinant associated with the autonomous system (17) is The advantage of the above equation is that it does not depend on the phase shift. Thus, for each pair of values k, η, (18) can be compactly rewritten as Solving (17) for e j φ in terms of V 1 and V 2 , one obtains And setting the magnitude of the expression in (20) to 1, one obtains the following system of three equations in three unknowns V 1 , V 2 , ω: where the superscripts r and i indicate real and imaginary parts. The practical solution of (21) is carried out in two stages. In a first stage, for each pair of values η, ω, one obtains two surfaces Then, one should calculate the intersection of the surface S 1 with the plane of zero value det r = 0, and the intersection of the surface S 2 with the plane of zero value det i = 0. For each η, ω, these intersections provide two curves C r (V 1 , V 2 ) and C i (V 1 , V 2 ) in the plane defined by V 1 and V 2 . Then, the potential solution points of (19) at the particular values η, ω are given by the intersections between the two curves C r and C i where P stands for point. For each η and sweeping ω, P η,ω (V 1 , V 2 ) gives rise to one or more curves in the plane defined by V 1 and V 2 , denoted by: D η (V 1 , V 2 ), where D refers to "determinant." However, equation (21)(c) must also be fulfilled.
As an example, Fig. 9(a) presents the curves D η (V 1 , V 2 ), composed by the points at which det(V 1 , V 2 , ω) = 0 obtained Fig. 9. Numerical method applied to the coupled transistor-based system in Fig. 1(b). Procedure to calculate the solution curves using conditions (18) and (24). for the coupling factor k = 0.21 and the capacitor value C 2 = 21 pF, which have been traced in blue. Note that all the points in the blue curves of Fig. 9(a) fulfill both det r = 0 and det i = 0. As stated, the actual solution points should also satisfy H = 1 [in (21)(c)]. For convenience, two functions are considered, respectively, obtained by solving for e j φ from each of the two complex equation composing the matrix system (17). These independent solutions are Note that when (18) is fulfilled, the two equations are identical. For each pair of values η, ω, by setting the magnitude of each of the two expressions in (23) to 1, one obtains one or more curves in the plane V 1 , V 2 . These curves are defined by When sweeping ω, the intersections of H 1 (V 1 , V 2 ) and H 2 (V 1 , V 2 ) provide one or more curves in the plane defined by V 1 and V 2 . These curves are denoted by: H η (V 1 , V 2 ). They have been traced in red Fig. 9(a). Solution points must correspond to intersections of D η and H η , which can be calculated thanks to the fine tuning enabled by the frequency sweep. In the plane defined by V 1 and V 2 , each point of the curves D η (V 1 , V 2 ) and H η (V 1 , V 2 ) corresponds to a distinct frequency ω. Only intersections at the same ω constitute solution points, which is verified through a simple error condition.
In the case of Fig. 9(a), only four intersections occur in the plane V 1 and V 2 between D η (V 1 , V 2 ) and H η (V 1 , V 2 ). Two of them correspond to actual solution points, as shown in Fig. 9(b) and (c), where D η and H η are traced in terms of V 1 and V 2 versus the frequency ω. The solution points are clearly distinguished. Now variations in the parameter η will be considered. In the case of the two coupled differential oscillators of Fig. 1(b), this parameter corresponds to the capacitance C 1 in the first oscillator. Fig. 10  The solution predicted by the numerical technique has been compared with experimental measurements. With this aim, the coupling factor k of the two coupled differential oscillators has been estimated from the measurement of the scattering matrix of the coupled inductors. The resulting value is k = 0.21. Fig. 11(a) and (b) presents a comparison of the predicted and measured variations of the oscillation amplitudes and phase shift versus the capacitor C 1 . The measurement points exhibit a good agreement with the analysis results corresponding to one of the oscillation modes. Fig. 11(c) presents the phase-noise spectrum obtained with the conversion-matrix approach [49] at C 1 = 22 pF. The results are compared with the experimental characterization using the R&S FSWP8 Phase Noise Analyzer. Both in simulations and measurements, the phase-noise spectrum is compared with the one obtained in uncoupled conditions (stand-alone operation). An improvement of about 3 dB is obtained, in agreement with the theory. The two periodic modes obtained for k = 0.21 are stable, with different basins of attraction. Fig. 12 presents timedomain simulations and experimental measurements, demonstrating the physical coexistence of the two modes. The basin of attraction of one of the modes is much larger than that of the others, so the out-of-phase mode was rarely observed in practice. Note that there can be situations in which only one of the modes is stable. This will happen, for instance, if only one the four sections (delimited by the turning points) of the original eight-shaped curve is stable.
V. COUPLING AT THE FREQUENCY RATIO 1:N The nonlinear analysis of two mutually injection-locked oscillators at the ratio 1: N is more involved than the one at 1:1, considered in Section III. To formulate the system, coupling effects are considered at the resonance frequency Nω, corresponding to the fundamental frequency of the higher frequency oscillator. Thus, the fundamental component of the first oscillator is only affected by the second oscillator through the coupling of its harmonic component V N at Nω. Assuming injection-locked operation at ω, the system is described as where Y T,1 is the total admittance function of the lower frequency oscillator at its fundamental frequency, Y A,N is the nonlinear active admittance function of this oscillator at the harmonic frequency Nω, Y A,b is the nonlinear active admittance function of the higher frequency oscillator at Nω, and y pi j are parameters of the admittance matrix of the coupling network, which may include additional linear elements. Its dependence on (ω, k, η) has been dropped for simplicity.  The phase origin is set at the Nth harmonic component of the first oscillator, so V N is real. The equation Y T,1 = 0 in (25)(a) can be taken as a constraint when addressing the coupled system (25)(b). This can be done by extracting the admittance function Y A,N (V 1 , V N , φ, ω) Fig. 14. Analysis of the system of two coupled differential oscillators at the frequency ratio 1:3 shown in Fig. 13. (a) Near conic surface in the space defined by V 1 , V 3 , and ω. (b) "Unfolded" surface in the space defined by V 3 , ω, and φ. under the fulfillment of Y T,1 = 0. If the oscillator at ω has good convergence properties in stand-alone (noncoupled) operation, the function Y A,N (V 1 , V N , φ, ω) can be calculated by introducing two AGs into this oscillator, one at ω (AG 1 ) and the other at Nω (AG 2 ). The AG 2 is not optimized but used as an excitation source to calculate the admittance  V N , φ, ω), calculated under the constraint Y T,1 = 0, will be introduced in (25)(b). If the oscillator at ω does not exhibit good convergence properties, the function Y A,N (V 1 , V N , φ, ω) can be obtained using the method in [26] for the calculation of synchronized-solution curves in a single injection-locked oscillator. Thus, extracting Y A,N (V 1 , V N , φ, ω) should not be a problem.
The procedure to solve (25)(b) is identical to the one described in Section III. The search for solution points of (25)(b) must be performed using the function Y A,N (V 1 , V N , φ, ω) calculated in the previous analysis and the pairs of values V N and ω must be limited to those resulting from that analysis. The values of V N and ω constitute a near conic surface. Instead, the points V b , ω have not limitation and define a plane.
The method has been applied to analyze the system in Fig. 13, coupled at the ratio 1:3. When isolated from each other, the two circuits oscillate at the respective frequencies 10 and 30 MHz. Fig. 14(a) presents the near conic surface in the space defined by V 1 , V 3 , and ω. The "unfolded" surface in the space defined by V 3 , ω, and φ is shown in Fig. 14(b). The solution curves obtained through (25) are shown in Fig. 15. In Fig. 15, the voltage amplitude in one of the oscillators, the synchronized oscillation frequency ω, and the phase shift φ have been represented versus the capacitance in the higher frequency oscillator for different k values. For the implementation shown in Fig. 13(b), the coupling factor estimated from the measurement of the scattering parameters of the coupled inductors is k = 0.15. Fig. 16(a) and (b) presents before and after synchronization when varying C 1 . Fig. 16(c) presents the phase-noise spectrum of the oscillator at 30 MHz, experimentally characterized with the R&S FSWP8 Phase Noise Analyzer, in both coupled and free-running operation.

VI. CONCLUSION
An in-depth investigation of the operation of two inductively coupled oscillators under strong coupling conditions has been presented. This is based on an analytical formulation of the coupled system, and the behavior pattern obtained when increasing the coupling factor has been found to be general and thus also observed in realistic transistor-based oscillators. Under the bilateral injection locking, one obtains two distinct families of solutions curves, one corresponding to each major operation mode. For the lower values of the coupling factor, one obtains a closed solution curve and two open curves, which for a higher coupling factor merge into two distinct curve families. A numerical method has also been developed for the realistic analysis of coupled transistorbased oscillators, which enables an exhaustive detection of all the coexisting oscillation modes. The method tackles the system equations in an original way that enables the use of contour-intersection techniques. The numerical method has been successfully extended to the case of oscillator coupling at the frequency ratio 1:N. This requires the extraction of the model of the lower frequency oscillation under the constraint of the fulfillment of the oscillation condition at its fundamental frequency. The methods have been applied to two pairs of coupled transistor-based oscillators at different frequency ratios and very good results have been obtained in comparison with the measurement results.