Chapter 8
Numerical Methods for Simulating Coupled-Oscillator Arrays
Coupled-oscillator arrays present a challenge to the designer due to difficulties both in the accurate simulation of oscillator elements and in the requirement for computationally efficient simulation techniques for large arrays. In addition, coupled-oscillator array design is made more difficult by the presence of multiple operating modes and stability considerations. As a result, a number of approximations need to be used to reduce the simulation time. These are describing function models for nonlinear elements [15, 118, 129, 130], along with perturbation models, infinite array approximations and continuum models shown in Chapter 3 [38, 39] and in Chapter 4 [42–44], respectively.
The progress in recent years in nonlinear simulation techniques has led to more accurate analysis and optimization methods for nonlinear circuits such as oscillators and mixers, as well as arrays [120,131]. Furthermore, these nonlinear simulation tools can be combined with electromagnetic simulation in order to analyze radiating structures and nonlinear antennas and arrays.
In this chapter, an introduction to numerical methods for simulating nonlinear circuits is presented [131, 132], focused on the simulation of autonomous circuits such as oscillators, followed by an introduction to convex optimization principles [133]. Nonlinear simulation techniques are demonstrated in order to trace the steady-state solutions of coupled-oscillator arrays and investigate their stability [116].
8.1 Introduction to Numerical Methods
The recent advances in numerical methods for simulating nonlinear microwave circuits permit one to model oscillator and coupled-oscillator array circuits efficiently and accurately. In this section a brief introduction to the principles of commonly used methods will be presented, with an aim toward obtaining the periodic steady state of oscillator circuits. The reader is prompted to the literature for an advanced and detailed description of the various methods, such as for example Refs. 101, 131, 132, and 134. Among the various existing numerical methods, transient simulation, harmonic balance, and envelope-transient simulation are described next.
8.1.1 Transient Simulation
A general nonlinear circuit is considered where a vector x of size N contains the state variables of the circuit, namely node voltages and currents. The circuit is described by a nonautonomous system of differential equations obtained by applying Kirchhoff’s current law at the circuit nodes as well as the voltage law at the circuit branches, as introduced in Eq. (7.1-2) and repeated here for convenience:
with an initial condition x(t0) = x0, where typically t0 = 0. A vector u of size P including external, known, forcing terms has been included for generality.
The system is classified as an initial value problem [132], and the computation of its solution over a given time interval is known as transient simulation. There exist various discrete time numerical integration methods that are used to perform a transient simulation [101]. Assuming an integration time step q, the values of the state variable vector xk at time tk = kq are generally computed as follows [101]:
The number of evaluations of the state variables and vector field that are required for the evaluation of the next state k + 1 are called the steps of the algorithm, and these steps define the order of the algorithm. An algorithm is called explicit if the future state depends only on past values of the state variables and the vector field, which corresponds to b−1 = 0. If b−1 ≠ 0, the algorithm is called implicit.
A commonly used single-step, explicit integration algorithm is the forward Euler algorithm, which is defined as
In contrast, the backward Euler algorithm is a single-step implicit algorithm
where the evaluation of the state vector at time k + 1 requires the computation of the vector field at the same time step. The resulting nonlinear system of algebraic equations maybe solved using some numerical root-finding algorithm. Typically the Newton–Raphson algorithm is used to compute the solution at each time step [132]. Assuming a nonlinear system h(x) of algebraic equations with unknown the steady state x = xk+1 at time step tk = (k + 1)q,
the Newton–Raphson algorithm is an iterative algorithm that requires an initial guess x0 as a starting point, and it proceeds to find the roots of h(x) by calculating successive approximations of the unknown steady-state vector as
where j is the iteration index, and jh(x(j)) is the Jacobian of the nonlinear function h(x(j)) [132]. The steady-state vector at the previous time step j is a good candidate for an initial guess x0 = xk. It can be shown that if the initial guess is close enough to a solution given by Eq. (8.1-5), if the nonlinear function h is continuously differentiable, and the Jacobian jh is not singular, the sequence given by Eq. (8.1-6) converges to a root of h.
Many different numerical integration algorithms (8.1-2) exist, depending on the choice of the various an and bn coefficients. Selection of the appropriate integration algorithm depends on computational complexity, accuracy, and numerical stability considerations [101]. Moreover, modern integration routines adaptively adjust the integration step and order of the integration algorithm.
In order to obtain the periodic steady state of an oscillator, one needs to integrate Eq. (8.1-1) for a sufficient time interval in order to allow all transient responses to decay. As a result, transient simulation is not an efficient method to analyze the behavior of oscillator and coupled oscillator systems. Conversely, however, transient simulation provides a way to examine the stability of the solutions, as time-domain integration converges only to stable steady-state solutions.
There exist time-domain algorithms such as the shooting methods that minimize the evaluation of the initial transient state in order to efficiently obtain the desired periodic steady state. In this case, one solves the system of differential equations given by Eq. (8.1-1) subject to a periodic boundary condition x(t) = x(t + T), where T is the period of the steady state. The reader is prompted to the literature for a description of these methods [101,132].
8.1.2 Harmonic Balance Simulation
Frequency-domain methods are particularly suited for the analysis of systems where a periodic solution exists. In this case it is possible to represent the steady-state solution by a trigonometric polynomial of degree M. The selected value of M is a trade-off between accuracy and computational efficiency. Specifically, assuming a state vector x(t) of size N, along with a vector of external forcing signals u(t) of size P, we can write
where ω is the angular fundamental frequency of the periodic steady state. The frequency-domain state vector and external forcing signal vector are defined by the N by 2M + 1 matrix X = [Xk] and by the P by 2M + 1 matrix U = [Uk], respectively. Similarly, the vector field f(x, u, t) is a periodic function and can also be expanded in a Fourier series as
where F = [Fk(X, U)] is the frequency-domain vector field N by 2M + 1 matrix, and it depends both on X and U.
In a typical piecewise harmonic balance algorithm implementation [134], the circuit is divided into a linear subcircuit and a nonlinear subcircuit, and Kirchhoff’s laws are applied in the nodes that connect the two subcircuits. The response of the nonlinear subcircuit is computed in the time domain, and a fast Fourier transform algorithm is used to convert the related data to the frequency domain. As a result, in order to compute the frequency-domain vector field matrix F, one first applies the inverse Fourier transform to the state X and external signal U vectors in order to obtain their time-domain expressions x and u, then computes the time-domain vector field f(x, u, t), and finally applies the Fourier transform to f(x, u, t) in order to obtain F.
By introducing Eqs. (8.1-7), (8.1-8), and (8.1-9) into the original time-domain system of differential equations given by Eq. (8.1-1) and then balancing the coefficients of the exponential terms ejkωt, a system of algebraic equations is obtained:
The matrix Ω contains the angular frequency terms generated by the time derivative operation on the Fourier series expansion in Eq. (8.1-7). The above system of algebraic equations is efficiently solved using root finding algorithms such as for example the Newton–Raphson algorithm [134,132] described in the previous section.
It should be noted that when the steady state is expanded using only a first-order trigonometric polynomial (M = 1), the corresponding formulation is known as the describing function [131], and it can be used to obtain insightful analytical expressions. The coupled-oscillator models of the previous chapter are describing function formulations.
Harmonic balance is able to handle quasi-periodic solutions by properly extending the polynomial basis and the time to a frequency-domain transform algorithm [132].
8.1.3 Conversion Matrix
In microwave mixer circuits, a quasi-periodic steady-state solution exists with two or more fundamental frequency components. In the simplest scenario, two fundamental frequencies need to be considered corresponding to the local oscillator signal and the RF input signal to the mixer. Correspondingly, a two-fundamental-frequency harmonic balance algorithm needs to be used in order to evaluate the steady state.
However, in typical mixer operation, the local oscillator signal has significantly larger power than the RF input to the mixer. As a result, it is possible to evaluate the periodic steady state in the absence of the RF input signal, defined by the local oscillator signal and using a harmonic balance algorithm with a single fundamental-frequency component. The effect of the RF input signal is then considered as a linear perturbation of the previously defined steady state, leading to a computationally efficient algorithm known as the conversion matrix method.
Assuming a linear perturbation of the steady-state solution x(t) = x0(t) + δx(t), along with an external RF signal u(t), the initial system of differential equations becomes
where x0(t) is the solution that corresponds to the large local-oscillator signal in the absence of the RF input, and dxf and duf are the time-varying Jacobians of the vector field f versus the state vector x(t) and versus the input RF signal vector u(t), respectively, evaluated at x0(t) and u(t) = 0. Both equations of the above system are solved in the frequency domain by applying the harmonic balance algorithm as Eq. (8.1-10). The frequency domain coefficients of the Jacobian matrices involved in the second equation are obtained at no additional computational cost during the Newton–Raphson harmonic balance computation of the large signal steady state corresponding to the first equation of Eq. (8.1-11) [131, 134].
8.1.4 Envelope Transient Simulation
The envelope transient simulation is a combination of the transient and harmonic balance simulation methods proposed by Sharrit [135] and by Ngoya and Larcheveque [136]. In effect, one represents the state variables, external forcing terms, and vector field by Fourier-series expansions of time-varying phasors:
Consequently, a transformed system of differential equations is obtained that has the form
The above system is solved using time-domain integration. The advantage of envelope transient simulation over the traditional transient simulation is that the time-varying phasors X(t) are slowly varying, allowing one to use a much larger time step in the simulation. Being a time-domain simulation, envelope transient simulation may also be used to verify the stability of a steady-state solution as it converges only to stable solutions.
8.1.5 Continuation Methods
Once a steady-state solution is obtained, continuation methods can be used to obtain the families of steady-state solutions that occur as one or more parameters of the circuit under consideration are varied. Continuation techniques provide an initial condition that is close to the required steady-state solution, so that the application of the Newton–Raphson or any other root finding algorithm that is being used converges quickly and efficiently.
Assuming a parameter p = p0 for which the steady-state solution x = x0 has been evaluated, it is then possible to obtain the steady-state solution xe corresponding to the parameter value p = pe by considering a sequence of values p0 < p1 < p2 … < pe and progressively evaluating the steady state corresponding to each parameter value by using the solution at each step as the initial condition for the evaluation of the next step [132].
In order to reduce the steps of the continuation method, the already obtained steady-state values are extrapolated. Assuming that the steady-state solution at step k has been obtained by solving the harmonic balance system [Eq. (8.1-10)]
the solution corresponding to pk+1 is approximated by linear extrapolation as
where JXHk is the Jacobian matrix of the harmonic balance system. The above matrix equation can be solved in order to obtain an initial condition for the state vector Xk+1:
Continuation methods based on Eq. (8.1-18) may fail due to singularities in the Jacobian matrix 
which result from the existence of multiple solutions versus the parameter under consideration. In this case, tracing of the steady-state solutions can be accomplished by parameter switching [137]. Parameter switching corresponds to tracing the steady-state solutions versus another, different circuit parameter or steady-state variable, for which the corresponding Jacobian matrix is not singular.
8.2 Obtaining Periodic Steady-State Solutions of Autonomous Circuits in Harmonic-Balance Simulators
Autonomous circuits, such as free-running oscillators, present an additional difficulty in harmonic-balance simulators due to the fact that the frequency basis of the trigonometric polynomial expansion is unknown. The autonomous nature of the oscillators is expressed in the time domain by the lack of a time reference, which translates in the frequency domain to an arbitrary phase of one of the harmonic components of its state variables. This fact is explored by Rizzoli and Neri [134], where it is proposed that the harmonic-balance system of algebraic equations can be extended by one more equation defined by arbitrarily setting the phase of one of the harmonics of a circuit state variable to a specific value. As an example, the phase of the fundamental-harmonic component may be set to zero leading to
Augmenting the harmonic-balance system by one equation allows one to additionally augment the number of unknowns by the fundamental frequency ω. Nonetheless, the Newton–Raphson algorithm may still converge to a DC (nonoscillating) solution due to the difficulty in selecting a suitable initial condition that is sufficiently close to the desired oscillating steady state.
Ch.-R Chang proposed an alternative method [138] in which an oscillator circuit is represented by a one-port equivalent circuit by looking into the terminals of the oscillator load, as shown in Fig. 8-1. The steady-state oscillation condition in the frequency domain is expressed by the total resistance or admittance at the load being equal to zero, known as the Kurokawa condition [129]. This condition, expressed at each harmonic k, is written as
Fig. 8-1. Oscillator 1-port equivalent circuit.
Enforcing the above condition in addition to the harmonic balance system enables the algorithm to avoid converging to solutions corresponding to zero harmonic components Vk and Ik, such as the nonoscillating DC solution. In fact, in order to avoid the DC solution, it is necessary to impose the admittance condition only at the fundamental harmonic component [138]
which leads to two additional real equations in the harmonic balance system
As a result, two additional variables can be introduced to the extended harmonic balance system: the unknown frequency ω and the oscillation amplitude V1 at the load. The additional advantage of this formulation is that the designer may impose in a circuit optimization problem the desired oscillation frequency and amplitude at the load. A dual formulation may also be obtained by considering a series one-port equivalent circuit and enforcing the oscillation condition by setting the total impedance equal to zero.
The condition given by Eq. (8.2-3) was implemented by Quere et al. in commercial simulators [139], allowing for a practical design and optimization methodology for autonomous circuits. According to Ref. 139, one needs to introduce into the simulator an ideal probe circuit such as the one shown in Fig. 8-2. The probe is connected in parallel to a selected circuit node and consists of an ideal sinusoidal source of a given amplitude Vs, phase 
s, and frequency fs, connected in series with a current meter Is and an ideal filter. The filter is such that it presents infinite impedance for frequencies other than the ideal source frequency fs, thus restricting the effect of the probe to fs.
In free-running oscillator simulation, the phase of the probe is set to an arbitrary but fixed value, for example zero. An optimization loop is run in order to find the nonzero amplitude and frequency of the probe that correspond to zero admittance Ys = Is/Vs. Each iteration of the loop is a harmonic-balance analysis. The result (Vs, fs) of the optimization defines the oscillating steady state. Alternatively, in the case of an externally injection-locked oscillator, the frequency fs is known and corresponds to the frequency of the external source. In this case, the pair (Vs, s) represents the unknowns of the optimization loop, as the oscillation phase is not arbitrary any more; rather, it depends on the injection source.
Fig. 8-2. Ideal probe circuit used for oscillator simulation in harmonic balance.
The ideal probe can also be used to initialize an envelope transient simulation to the oscillating steady state [140]. The optimization loop is first run in order to obtain the oscillating steady state (Vs, fs), and subsequently an envelope-transient analysis is executed with the probe connected to the circuit only for an initial small time interval. The probe is then disconnected from the circuit (for example with the help of a time-dependent switching resistor), and the circuit is left to evolve for the remaining time interval according to its dynamics. This way, the envelope-transient analysis can be used to verify the stability of the steady-state solution. Once the probe is disconnected from the circuit, if the solution is unstable, the circuit will evolve to a different steady state.
8.3 Numerical Analysis of a Voltage-Controlled Oscillator
The simulation tools described in the previous section are now used to design a voltage-controlled oscillator that may serve as the array element in the coupled-oscillator array numerical analysis examples of the following sections. The oscillator circuit is based on the pseudomorphic high electron mobility transistor (pHEMT) device shown in Fig. 8-3.
A series resonator is connected at the gate terminal of the device, and a feedback capacitor is introduced at the source terminal. The feedback capacitance guarantees the presence of a negative resistance at the gate terminal. At the output, a matching network composed of two inductors is formed at the drain terminal. A frequency-tuning varactor is connected at the source terminal. The pHEMT device is self-biased, using a resistor placed at the source terminal of the device. Additionally, a 50-ohm (Ω) termination is used at the gate terminal in order to accommodate a port for an external injection signal, and also bias the gate terminal at 0 V DC, ensuring a negative gate-source (VGS) voltage. The values of the various circuit components and bias are shown in Table 8-1.
Fig. 8-3. Oscillator circuit schematic.
Table 8-1. Oscillator circuit component values.
| Parameter | Value |
| Ldd (nH) | 0.5 |
| Ld (nH) | 0.2 |
| Cd (pF) | 1.0 |
| Lg (nH) | 3.3 |
| Cg (pF) | 0.5 |
| Cs (pF) | 1.5 |
| Rc (k Ω | 5.0 |
| Ri (Ω) | 50.0 |
| RL (Ω) | 50.0 |
| Rs (Ω) | 25.0 |
| Vdd (V) | 1.5 |
| Vc (V) | 0-15 |
| pHEMT | NE3210S01 |
| D1 | MA46H070-1056 |
nH = nanohenry, pF = picofarad
Harmonic balance optimization using an ideal probe to ensure convergence to the oscillating steady state was used for the design. The probe was connected to the output node ν0; however, other nodes may also be used such as any of the pHEMT terminals. The use of the output node is convenient because one can directly optimize the oscillator output power.
The VCO frequency and output power are shown in Fig. 8-4 and Fig. 8-5, respectively. The oscillator is consuming about 22.5 mW from a 1.5-V supply for all values of the control voltage. Its DC-to-RF conversion efficiency is approximately 9.5 percent at Vc = 0 V and reduces to 6.3 percent as Vc reaches 15 V.
The ideal probe is then used to simulate oscillator synchronization properties when an external injection source is applied at the gate terminal. In the first case, an external source with fixed power is connected to the circuit and its frequency is varied in order to obtain the synchronization curves. The control voltage is fixed to 10 V. The result is shown in Fig. 8-6 for two values of available power of the injection source. The synchronization bandwidth is equal to the frequency interval contained between the two edges of the closed curves defined by the infinite slope of the power versus frequency curves. The free-running steady state is represented by a point in the plot corresponding to a frequency of 9.892 GHz and power of 2.9 decibels referenced to milliwatts (dBm). As the injection power increases, the synchronization curves become larger, and they eventually open [141].
Fig. 8-4. VCO frequency versus the control voltage.
Fig. 8-5. VCO output power versus the control voltage.
Fig. 8-6. VCO synchronization curves versus the injection signal frequency for a fixed control voltage Vc = 10 V. The free-running frequency and power are indicated by a point ‘x’ in the plot.
Alternatively, one may fix the injection signal frequency and obtain the synchronization curves versus the control voltage Vc, which corresponds to varying the free-running frequency of the oscillator. The corresponding curves obtained for an injection-signal frequency of 9.892 GHz (equal to the free-running oscillator frequency at Vc = 10 V) are shown in Fig. 8-7. It should be noted that the curves of Fig. 8-6 and Fig. 8-7 are generated by sweeping the phase difference between the injection signal and the oscillator output by 360 deg.
The points of infinite slope are the turning points of the curve, and due to the fact that the synchronization curves are closed, for any frequency or control voltage between them, there exist two solutions for the oscillator power.
In Fig. 8-7, the free-running frequency and power are indicated by a point ‘×’ in the plot. For every Vc value between the turning points, two steady-state solutions exist (for example points A1 and A2 correspond to Vc = 10.45 V).
It was shown in Section 7.2 that turning points correspond to a change of the stability of the steady-state solution; as a result, only one of the two solution branches joined by the two turning points is stable. Each solution branch corresponds to a phase shift variation of 180 deg between the injection signal and the oscillator output. Therefore, in the case of an injection signal with frequency close to the fundamental frequency of the oscillator, one can obtain approximately up to 180 deg of (stable) phase shift range between the injection signal and the oscillator signal.
Fig. 8-7. VCO synchronization curves versus the control voltage for a fixed injection-signal frequency finj = 9.892 GHz, equal to the oscillator free-running frequency at Vc = 10 V.
It is possible to determine which one of the branches is stable by selecting one point on it and running an envelope transient analysis initialized to the steady-state solution as was described in Section 8.2. For example, for an injection signal power of −35 dBm, the points A1 and A2 of output power 2.35 dBm and 3.4 dBm, respectively, correspond to Vc = 10.45 V (Fig. 8-7). The stability of A2 was examined by envelope transient simulation. A simulation time of 100 nanoseconds (ns) and a time step of 5 ns were used, with a Gear time-domain integration algorithm [101]. The result is shown in Fig. 8-8, where one can see that the system evolves to point A1, demonstrating that the upper branch of Fig. 8-7 is stable. In a similar way, one can also verify that the upper branch of Fig. 8-6 is the stable one.
Finally, it should be noted that the oscillator admittance derivatives of the perturbation model of Section 7.6 that is used to model coupled oscillator arrays, can be easily computed from a harmonic-balance simulation of the single element with an ideal probe placed at the desired oscillator node. The oscillator circuit admittance derivatives are equal to the derivatives of the ideal probe admittance. Once the steady state corresponding to zero admittance looking into the probe has been determined, the probe admittance derivatives can be evaluated using finite differences [116]. As an example, the admittance derivatives for the oscillator of Fig. 8-3, corresponding to the steady state defined by control voltage Vc = 10 V, frequency f0 = 9.892 GHz, and amplitude 0.442 V (P0 = 2.9 dBm), are listed in Table 8-2.
Fig. 8-8. Envelope-transient analysis of the steady state corresponding to A2 of Fig. 8-7. A2 is unstable and the system evolves to A1.
Table 8-2. Oscillator steady-state and corresponding admittance partial derivatives calculated using finite differences
| Parameter | Value |
| P0 (dBm) | 2.9 |
| F0 (GHZ) | 9.892 |
| Vc (V) | 10.0 |
 |
0.0547 + j 0.1957 |
 |
0.002 − j 0.008 |
 |
-0.015 + j 0.008 |
8.4 Numerical Analysis of a Five-Element Linear Coupled-Oscillator Array
The VCO of the previous section is used here to create a five-element linear coupled-oscillator array. The array elements are coupled with resistor-loaded transmission-line sections of 50-Ω characteristic impedance and electrical length of 360 deg at a frequency of 9.89 GHz, which corresponds to a control voltage of 10 V in the free-running VCO element. The series resistors in the transmission line coupling sections control the coupling strength among the array elements, as was proposed by Liao and York [142]. The schematic of the array is shown in Fig. 8-9. The coupling network is connected at the oscillator outputs and each oscillator can be externally injected through its gate terminals.
Fig. 8-9. Five-element linear coupled oscillator array.
Harmonic-balance simulation is used to trace the various solutions of the coupled oscillator array, by connecting one oscillator probe at the output node of each oscillator element. That way, it is guaranteed that the simulator will properly converge to the periodic steady state of each oscillator in the array. The five probes extend the harmonic balance system of algebraic equations by ten real equations, thereby allowing the designer to optimize ten additional unknowns. The synchronized solutions (which correspond to a constant phase shift among adjacent oscillator elements) are obtained by sweeping the phase shift while optimizing the five oscillator-output voltages, the common oscillation frequency, and four control voltages, all except the one corresponding to the middle oscillator.
The simulation results are shown in Figs. 8-10 through 8-12, where the output power, the frequency, and the control voltages, respectively, are plotted versus the oscillator phase shift. The coupling-network resistor is set to R = 270 Ω, and the control voltage of the middle oscillator is fixed at Vc3 = 10 V. The phase shift has been swept from 0 to 180 deg, with the oscillator phases increasing from oscillator 1 to the left and toward oscillator 5 to the right of Fig. 8-9. Due to the symmetry of the array, the solution curves for the remaining phase-shift values (0 to −180 deg) can be obtained by considering the mirror image of the array elements with respect to the central element 3—in other words, replacing element 5 with element 1 and replacing element 4 with element 2.
Figure 8-12 shows the variation of the oscillator-control voltages versus the phase shift. One can see that the edge element-control voltages present a significantly larger variation compared to the inner elements. In fact, the control voltages of elements 2 and 4 remain practically constant for phase shifts up to 90 deg. This represents a numerical verification using a harmonic balance simulation of the proposition of Liao and York [142], where, by only tuning the free-running frequency of the peripheral elements of a coupled-oscillator array, it is possible to generate constant phase-shift distributions among the array elements, thus both minimizing the required number of controls and eliminating the need for phase shifters.
Fig. 8-10. Five-element linear coupled-oscillator array. Output power of each oscillator versus the phase shift among adjacent elements. The coupling resistor is R = 270 Ω.
Fig. 8-11. Five-element linear coupled-oscillator array. Array frequency versus the phase shift among adjacent elements. The coupling resistor is R = 270 Ω.
Fig. 8-12. Five-element linear coupled-oscillator array. Control voltage of the oscillator elements versus the phase shift among adjacent elements. The middle oscillator element control voltage is fixed at Vc3 = 10 V, and the coupling resistor is R = 270 Ω.
However, one can observe from Fig. 8-11 that as the phase shift varies, the array frequency also varies. One way to maintain a constant frequency for every desired phase shift is by allowing the control voltage of one more oscillating element to vary. In this example, one may allow the middle oscillator-element control voltage (Vc3) to vary, thus being able to eliminate potentially undesired frequency variations. In Figs. 8-13 and 8-14, the five oscillator amplitudes and control voltages are plotted versus the phase shift for a coupling resistor of R = 330 Ω, while the array frequency is fixed at 9.892 GHz. The inner-oscillator control voltages take very similar values; however, they need to be varied in order to maintain the frequency of the array constant.
The coupling strength among the oscillator elements is set by the coupling resistor R. In fact the selection of the optimum coupling strength is a trade-off among a number of parameters. As the coupling strength increases, the variation in the output power of the oscillators, the frequency, and the control voltages with the phase shift all increase. As the coupling strength decreases, the oscillators eventually desynchronize due to the presence of noise.
Fig. 8-13. Five-element linear coupled-oscillator array. Output power versus phase shift for a coupling resistor R = 330 Ω. The array frequency is fixed at f = 9.892 GHz, and the control voltages of all elements are allowed to vary.
Fig. 8-14. Five-element linear coupled-oscillator array. Oscillator control voltages versus phase shift for a coupling resistor R = 330 Ω. The array frequency is fixed at f = 9.892 GHz.
It is easily verified from Figs. 8-10 to Fig. 8-12 that for the value of the coupling resistor of R = 270 Ω, a harmonic balance solution for every possible phase shift value exists. As the coupling strength increases, it is not possible to obtain a solution for every phase shift. This is demonstrated in Fig. 8-15, where oscillator output power of the harmonic balance solutions corresponding to a coupling resistor of R = 220 Ω is plotted. It is easily seen that solutions exist only up to approximately 120 deg, and they are limited by the presence of a turning point [116]. As the desired phase shift progressively increases, the amplitude of oscillator 2 eventually drops to zero.
The stability of the solutions was also verified using envelope transient analysis, using the method described in Section 8.2. The simulation results showed that the coupled oscillator array with R = 270 Ω loses stability for phase shift values larger than approximately 58 deg. Additionally, the perturbation model of Section 7.6 was used to evaluate the constant phase shift steady-state solutions and their stability. The real part of the largest eigenvalue of the linear variational equation of the steady-state solution is shown in Fig. 8-16, where one can see that the perturbation model predicts loss of stability for a phase shift approximately equal to 63 deg, a value that is in relatively good agreement with the result obtained from envelope transient simulation.
Fig. 8-15. Five-element linear coupled-oscillator array. Output power versus phase shift for a coupling resistor R = 220 Ω. As the coupling strength increases solutions do not exist for every phase shift value. The middle oscillator element control voltage is fixed at Vc3 = 10 V.
Fig. 8-16. Stability analysis of the steady-state solution using the perturbation model of Section 7.8 showing critical eigenvalue real part versus the phase shift between adjacent elements.
In Fig. 8-16, positive values of the real part correspond to unstable solutions. The coupling resistor is R = 270 Ω, and the control voltage of the middle oscillator is fixed at Vc3 = 10 V.
Further comparison between the harmonic balance solution and the perturbation model is made in Figs. 8-17 and 8-18, where the amplitude of oscillators 1 and 3 along with the array frequency are plotted versus the phase shift between adjacent elements obtained using both methods. One can see that the agreement becomes worse for large phase offsets where the perturbation is larger. The perturbation model is limited to small perturbations around the free-running steady state, which in this case is near the 0-deg phase shift (in-phase) solution, and, additionally, to oscillators with small harmonic content [116,143]. Nonetheless, the advantage of the perturbation model lies in its computational efficiency which quickly becomes important as the number of array elements increases.
Fig. 8-17. Comparison of harmonic balance simulation and perturbation model of Section 7.8. Output power of oscillators 1 and 3 versus the phase shift between adjacent elements. The coupling resistor is R = 270 Ω, and the control voltage of the middle oscillator is fixed at Vc3 = 10 V.
Fig. 8-18. Comparison of harmonic-balance simulation and perturbation model of Section 7.8. Array frequency versus the phase shift between adjacent elements. The coupling resistor is R = 270 Ω, and the control voltage of the middle oscillator is fixed at Vc3 = 10 V.
8.5 Numerical Analysis of an Externally Injection-Locked Five-Element Linear Coupled-Oscillator Array
Injection locking the array to an external source signal is desirable in several applications in order to reduce the array phase noise as shown by Chang et al. [123], or introduce modulation to the oscillator signal as considered by Kykkotis et al. [99] and Auckland et al. [122].
The dynamics of the system and the stability of the various solutions depend strongly on the element that is being injected, whether it is located in the center of the array or near the edges [120, 123, 140]. Additionally, the number of elements that are being externally injected strongly influences the number and behavior of the existing solutions. Commonly used topologies are the one proposed by Stephan [1], in which the two end elements of a linear array are injection-locked to an external source, and the topology where the external signal is illuminating all the elements of the coupled oscillator array leading to a globally injection locked array [123], such as the case of a reflectarray or transmit-array antenna.
In the case of an externally injection-locked array, the oscillation frequency is determined by the frequency of the external source. In contrast, the phase difference between the injection source and the element that is being injected must be included in the unknowns of the harmonic-balance system of equations. Similarly with the free-running array case, a probe must be connected to each oscillator element in order to guarantee the convergence of the harmonic balance simulator to the oscillating solution.
The five-element array of Section 8.5 is considered with a coupling resistor of R = 330 Ω. The middle element (3) is injection-locked to an external signal source through its gate termination. The steady-state solutions corresponding to a constant phase shift among the array elements are traced versus the phase shift among adjacent elements. The additional unknowns in the harmonic balance optimization that can be obtained due to the use of the ideal probes are the five oscillator amplitudes, the four control voltages corresponding to all the elements (except the one being injected), and the phase difference between the injected element and the external source signal. The phase of the injected element is fixed at 0 deg, and the phase of the injection signal inj is allowed to vary. The control voltage of the injected element is fixed at 10 V. Finally, the frequency of the external signal is 9.892 GHz.
In Fig. 8-19 the phase inj is plotted versus the phase shift between the oscillator elements for different injection-signal powers. As was the case in Section 8.4, due to the symmetry of the array, solution curves also exist for the phase shift interval between 0 deg and −180 deg, and they can be obtained by taking the mirror image of the array elements with respect to the central element 3. One can see that for a given injection-signal power and oscillator phase shift it is possible to have two solutions corresponding to two different values of inj.
The output power of the middle oscillator is plotted in Fig. 8-20 versus the inter-oscillator phase shift, for different injection signal-power levels. It can be seen that for lower injection-signal power levels (Pinj = −35 dBm), solutions for every inter-element phase shift in the range [0 deg, 180 deg) do not exist. Specifically, a closed solution curve exists for phase shifts up to approximately 100 deg, where a turning point appears. As the injection power increases, the solution curve widens, and the turning point shifts to a larger phase-shift value. Finally for larger injection powers, the curve opens, forming two solution branches. As shown in Fig. 8-20, as many as three solutions may exist for a given phase shift value.
Fig. 8-19. Externally injection-locked five-element linear coupled-oscillator array. Injection-signal phase versus the oscillator phase shift, for different injection-signal power levels. The injection signal frequency is 9.892 GHz, and the coupling resistor R = 330 Ω. The middle oscillator element-control voltage is fixed at Vc3 = 10 V.
Fig. 8-20. Externally injection-locked five-element linear coupled-oscillator array. Middle oscillator output power versus the phase shift, for different injection signal power levels. The injection signal frequency is 9.892 GHz, and the coupling resistor R = 330 Ω. The middle oscillator element control voltage is fixed at Vc3 = 10 V.
Subsequently, the stability of the solutions must be examined, in order to determine which of the multiple solutions are stable and will appear in practice. The solution stability may be determined using transient or envelope-transient simulation, or by examining the eigenvalues of the linear vanational system of equations corresponding to one of the analytical model formulations of the previous sections. In Ref. 120, Collado and Georgiadis studied the injection-locked solutions of a two-element array, and it was determined that there exists an optimum coupling strength that leads to a maximum stable constant phase-shift range.
8.6 Harmonic Radiation for Extended Scanning Range
When an array of oscillators is coupled at the fundamental frequency, the maximum stable phase-shift range that can be introduced between adjacent array elements is 180 deg. In the case of a coupling network that has a coupling phase of 0 deg, this translates to being able to generate constant phase-shift distributions Δ in the −90-deg to 90-deg range. Considering a linear array where the radiating elements are placed at a half-wavelength distance, its main beam can be scanned according to Δ = kd sin θ = π sin θ, where θ is measured from broadside, for a maximum of θ = 30 deg.
It is possible to extend the phase-scanning range by considering the fact that the phase variation of the oscillator Nth harmonic is N times the phase variation of its fundamental frequency component, where N is the harmonic order. The task of the designer then becomes that of being able to generate sufficient power in the desired harmonic component. Essentially there are two ways of implementing such architectures, either by placing a frequency N-tupler circuit at the output of each oscillator or by properly designing the oscillator elements to have maximum power at the harmonic under consideration.
In Ref. 145, Alexanian et al. proposed a linear array of five coupled oscillators, where each oscillator element is followed by a frequency doubler, as shown in Fig. 8-21. The fundamental frequency of the oscillators is 4 GHz, and their output power is 9 dBm. The prototype array in Ref. 145 used a compact field-effect transistor (FET)-based frequency doubler circuit with 1-dB conversion gain. The theoretical phase-tuning range that can be achieved with this topology is 360 deg.
Based on the same principle, a frequency-tripled two-dimensional coupled-oscillator array operating in X-band was reported by Pogorzelski in Ref. 69. An inter-oscillator phase difference ranging up to 60 deg was tripled to 180 deg. Thus, this array had a demonstrated H-plane scanning range of ±90 deg. The fabricated prototype additionally contained a diagnostic system used to evaluate the phase differences between the various oscillator elements. The array is described in more detail in Section 6.2.
Alternatively, Sanagi et al. [146] proposed a four-element coupled-oscillator array, where the oscillator elements were specifically designed in order to have a high second-harmonic content, thus also obtaining a 360-deg phase-scanning range. The proposed circuit is shown in Fig. 8-22.
The oscillators are coupled using directional couplers. Termination circuits based on the coupler networks are also attached to the edge elements in order to implement a symmetrical coupling network where all oscillators see approximately the same load. Sanagi et al. [146] extended the coupled oscillator model based on the cubic nonlinearity, which was introduced by York [111], in order to study their proposed circuit architecture. Specifically, they considered a nonzero square term in the cubic polynomial describing the current-to-voltage characteristic of the nonlinear device used for the oscillators; additionally, Sanagi et al. [146] introduced in the formulation an additional equation pertaining to the second harmonic. The block diagram of the considered circuit topology is shown in Fig. 8-23, which was used to investigate the effects in the array performance due to coupling both at the fundamental frequency and at the second harmonic. It was shown that as the second harmonic coupling becomes stronger relative to the coupling at the fundamental frequency, the achievable phase tuning range is reduced.
Fig. 8-21. Coupled-oscillator array using frequency doublers for extended scanning range. (Reprinted with permission from Ref. 145, ©1995 IEEE.)
Fig. 8-22. Coupled-oscillator array radiating the second harmonic frequency component (dimensions in millimeters). (Reprinted with permission from Ref. 146. This material is reproduced with permission of John Wiley & Sons, Inc.)
Fig. 8-23. Model of the coupled-oscillator array radiating the second harmonic-frequency component. (Reprinted with permission from Ref. 146. This material is reproduced with permission of John Wiley & Sons, Inc.)
In Ref. 147, Georgiadis proposed a three-element coupled-oscillator array shown in Fig. 8-24, also optimized in order to radiate the second-harmonic output wave. The array operates based on the same principle as the one by Sanagi et al. [146]. In this work, however, the alternative perturbation model for the coupled oscillator array given in Section 7.6 was extended in order to include the formulation for the second-harmonic frequency component.
Furthermore, harmonic-balance analysis was used to trace the steady-state solutions corresponding to constant phase shifts between the array elements at the second-harmonic component. In order to do so, it is necessary to place two ideal probes at each oscillator output node, one at each harmonic [147]. The computational load associated with the optimization of the coupled-oscillator array radiating the second harmonic is increased due to the fact that the number of ideal probes, and therefore optimization goals, required for the simulation is doubled.
8.7 Numerical Analysis of a Self-Oscillating Mixer
Self-oscillating mixers (SOMs) are particularly attractive for low-cost, compact implementations of microwave circuits due to the fact that the same circuit is used to provide a local-oscillator signal as well as for frequency translation.
Fig. 8-24. Three-element coupled-oscillator array prototype, designed to radiate the second harmonic frequency component. (Reprinted with permission from Ref. 147, ©2007 IEEE.)
The performance parameters of self-oscillating mixer circuits (such as conversion gain and inter-modulation distortion) can be evaluated using harmonic-balance simulation, provided that an ideal probe is used to enforce the convergence of the simulator to the oscillating steady state. The probe equations are set up in order to make sure that the admittance associated with the probe is equal to zero at the oscillating frequency of the circuit.
The radio-frequency (RF) and intermediate-frequency (IF) signals are treated by introducing a second fundamental frequency component in the harmonic-balance frequency basis, thus using a two-fundamental harmonic balance system of equations. Alternatively, one may consider the effect of the RF and IF signals as a linear perturbation of the oscillating steady state and employ the conversion matrix method to efficiently compute the conversion gain of the self-oscillating mixer. Finally, the RF and IF frequency signals can be efficiently treated using an envelope-transient simulation that has been initialized to the oscillating steady state.
In Ref. 148, Herran et al. optimized the gain associated with a selected mixing product of a self-oscillating mixer by using two ideal probes properly introduced in the circuit and optimizing the reflection coefficients of an ideal multi-harmonic load connected to the circuit input. The circuit schematic that was used is shown in Fig. 8-25.
The first probe, called an auxiliary generator in Fig. 8-25, is used to enforce the oscillation condition at the desired frequency. The admittance looking into this probe is set to zero in order not to perturb the circuit steady state, and the complex admittance or reflection coefficient of the multi-harmonic load at the fundamental frequency that satisfies this condition is found through harmonic balance optimization.
The second ideal generator probe is connected in series with the gate terminal of the FET device, and its frequency corresponds to a desired Nth harmonic that is selected for the mixing process. Mixing products involving the second and third harmonics were considered. The reflection coefficient of the multi-harmonic load at the desired harmonic is set to −1, corresponding to a short circuit. The optimization procedure consists of finding the complex amplitude of the ideal generator which results in a desired mixing gain value. The corresponding admittance looking into the generator must have a positive real part in order for it to correspond to a passive load. In this way, the multi-harmonic load is optimized for a desired mixing gain value and its reflection coefficient at the fundamental frequency and selected harmonic frequency are determined. The final design is obtained by implementing the obtained reflection coefficient values using passive printed or lumped circuit components [148].
A varactor diode may be appropriately placed in the self-oscillating mixer circuit in order to provide a frequency-tuning capability. An externally injection-locked self-oscillating mixer operates as both a mixer and a phase-shifter element, where the phase shift between the input and output of the mixer is varied by changing the free-running frequency of the self-oscillating mixer.
Being a synchronized oscillator, the externally injection locked self-oscillating mixer can be used to provide a continuous phase-shift range N × 180 deg, where the external injection signal is assumed to have a frequency near the fundamental frequency of oscillation of the self-oscillating mixer, and the Nth oscillator harmonic is used in the mixing operation. Here, the fact that the tuning range of the phase of the oscillator Nth harmonic is N times the tuning range of the phase of its fundamental frequency component being used [145].
Fig. 8-25. Nonlinear optimization of a self-oscillating mixer, (a) Circuit topology: The input signal of power Pin passes through a band-pass filter with center frequency fin, before it is mixed by the active circuit and collected at the output through an intermediate frequency filter of center frequency fIF. The optimization procedure consists of designing a multi-harmonic load with impedance Zn at harmonic nω0. Optimization is performed using an ideal auxiliary generator probe AG with amplitude VAG and frequency ωAG defined as in Fig. 8.2, along with a substitution generator with amplitude 
at the harmonic frequency Nω0. (b) Definition of the multi-harmonic load using ideal circuit components. The figure indicates the input reflection coefficients Γn corresponding to the load impedance Zn at frequency nω0. (c) Implementation of the multi-harmonic load using microstrip components for the cases of a second (2 HSOM) and third (3 HSOM) harmonic self-oscillating mixer, respectively. (Reprinted with permission from Ref. 148, ©2006 IEEE.)
It has been argued in Sections 7.9 and 8.5 dealing with the analysis of externally injection locked oscillator arrays that such architectures can be used to transmit information by introducing phase or frequency modulation in the external injection signal. Furthermore, the effect of modulation in the array scanning range was investigated in Section 7.11. Such topologies are limited to relatively narrowband applications due to the fact that the modulation strongly affects the steady state of the synchronized oscillator signals. Furthermore, specific modulation formats leading to small envelope variations are suitable for such applications due to the fact that the amplitude-limiting properties of the oscillators tend to introduce distortion to the envelope of the modulating signals. Continuous phase modulation (CPM) [149], which is a constant envelope modulation, is a prominent candidate for such systems. A well-known example of CPM is Gaussian minimum-shift keying (GMSK) used in the Global System for Mobile Communications (GSM), second-generation mobile (cellular) communication systems.
However, when modulation is introduced through the RF input signal of the self-oscillating mixer, it does not strongly affect the synchronization state of the mixer due to the fact that the input signal has a low power level and represents only a perturbation of the steady state. As a result, self-oscillating mixers can be used as frequency translation and phase-shifter circuits for input RF signals of arbitrary modulation. Furthermore, proper design of the mixer can allow one to obtain broadband gain, and therefore the self-oscillating mixer is not limited to RF input signals with narrowband modulation.
The use of an injection-locked self-oscillating mixer as a downconverter and phase shifter element was studied by Ver Hoeye [80]. The proposed circuit topology is the same as in Fig. 8-25 with the addition of a varactor diode connected in parallel with the series feedback shorted stub present at the source terminal of the active device in order to provide a frequency tuning capability. The SOM design was performed using the methodology described previously in this section. An oscillation at 3.25 GHz was obtained, and an RF signal of 11.25 GHz was mixed with the third harmonic of the SOM, resulting in an IF output of 1.5 GHz. Phase tuning of as much as 3 × 180 deg = 540 deg was achieved by utilizing the third harmonic mixing product. The obtained conversion gain was 4.5 dB over a bandwidth of approximately 100 MHz. It is shown in Fig. 8-26 that the conversion gain depends both on the injection power level Ps and on the varactor control voltage Vcont or, in other words, the selected phase difference between the input and output SOM terminals. The results have been obtained using a two-fundamental-harmonic balance simulation, and one can observe the closed synchronization curves of the injection locked self-oscillating mixer, which are similar to the ones obtained for the synchronized oscillator in Fig. 8-7. The synchronization curves open as the injection power increases, and there exist two solutions for a given control voltage within the synchronization band limited by the curve edges of infinite slope. Only one of the two solutions is stable and therefore measured experimentally, and in this case it corresponds to the branch with lower conversion gain.
8.8 Conclusion
In this chapter we provided an introduction to nonlinear analysis methods with a special focus on methods of analysis applied to nonlinear circuits such as oscillators, self-oscillating mixers, and coupled oscillator arrays. Such tools can be combined with electromagnetic simulators in order to accurately model the various passive components of the circuits under consideration such as transmission lines, interconnects, resonators, and antennas. Typically, these methods can be used to analyze small arrays consisting of tens of elements or fewer due to their increased computational complexity. Additionally, they can be used to compute the various parameters that are required to formulate the approximate models of the previous chapter such as the nonlinear admittance derivatives, which, in turn, can be used for an efficient less time-consuming simulation and optimization of the arrays.
Fig. 8-26. SOM conversion gain versus the varactor control voltage for different injection power levels Ps. (Reprinted with permission from Ref. 80, ©2006 IEEE.)