Understanding Relaxation Oscillator Circuits Using Fast-Slow System Representations

We demonstrate the utilization of the fundamental principle of nonlinear dynamics, namely, the Liénard-type representations of ordinary differential equations, also referred to as fast-slow systems, to describe and understand relaxation oscillations in electronic circuits. Relaxation oscillations are characterized by periods of slow signal changes followed by fast, sudden transitions. They are generated either intentionally by means of usually simple circuits or often occur unintentionally where they would not have been expected, such as in circuits with only one dominant energy storage device. The second energy storage required to promote oscillatory solutions of the governing equations can also be provided by spurious elements or mechanisms. The conditions that distinguish harmonic from (anharmonic) relaxation oscillations are discussed by considering the underlying eigenvalues of the system. Subsequently, we show how to intuitively understand relaxation oscillations through analyses of the phase diagram based on the fast-slow system representation of the nonlinear differential equation. Practical examples of oscillators including $RC$ and $LR$ op-amp circuits and the so-called “Joule thief” circuit are discussed to illustrate this principle. The applicability of the method is not limited to electrical circuits, but extends to a variety of disciplines, such as chemistry, biology, geology, meteorology, and social sciences.


I. INTRODUCTION
Relaxation oscillators are known to electrical engineers as simple circuits and are commonly taught in undergraduate electronics courses.Fig. 1 shows one of these circuits, which utilizes an op-amp configured as an inverting Schmitt trigger and an RC circuit feeding the output back to the input.The principle of operation is explained by considering the capacitor-charging equation which applies when capacitor C, with initial voltage V C,0 , is charged through resistor R. V C,∞ denotes the terminal voltage reached after infinite time.We assume that initially V C = 0, and the output voltage of the op-amp is saturated to the maximum positive output swing V B .The capacitor charges until V C reaches the upper switch threshold V T = The associate editor coordinating the review of this manuscript and approving it for publication was Derek Abbott .
V B R 1 /(R 1 +R 2 ).At this point, the output voltage swings to the negative rail −V B by virtue of the positive feedback and subsequently discharges the capacitor until the negative switch threshold −V T is reached.This, in turn, causes the output to swing back to the positive supply rail.The resulting output signal is a square wave with a 50 % duty cycle.The period T is obtained using (1), by setting V C (T /2) = −V T , V C,0 = V T and V C,∞ = −V B and solving for T yielding where log denotes natural logarithm.The term ''relaxation'' refers to the fact that the output signal shows only small, or in this case, no changes over a certain time duration but then suddenly performs a fast change, i.e., it ''relaxes''.While this behavior is not surprising for the shown circuit due to the hysteresis of the Schmitt trigger, relaxation oscillations also occur in circuits where no obvious discontinuous functions are involved; however, hysteresis is found as a result of some kind of nonlinear behavior.
(for a rail-to-rail op-amp).The capacitor voltage V oscillates between ±V T giving a square output signal V O .
The first formal investigation and the coining of the term ''relaxation oscillation'' can be traced back to the 1920s, when Dutch electrical engineer Balthasar van der Pol observed these phenomena in a triode vacuum tube circuit, as shown in Fig. 2 [1], [2], [3].Although the vacuum tube circuit was expected to exhibit harmonic oscillations owing to the soft nonlinear transfer characteristic of the vacuum tube and the linearity of the remaining circuitry, van der Pol noticed the presence of relaxations.His formal analysis led to a deeper understanding of these phenomena, paving the way for understanding relaxation oscillations in various fields, including mechanics, biology, chemistry, and engineering problems [4].These oscillations have been found in diverse systems such as nerve activity and heart beat [5], the Zhabotinsky reaction [6], earthquakes, prey-predator and business cycles [7], [8], climate-related cycles [9], stick-slip oscillations of bowed violin strings [8], and memristor oscillators [10], [11].It is fascinating that simple circuit theory and basic concepts of dynamical systems are sufficient to understand these far-reaching phenomena.As such, engineers familiar with these concepts hold the key to unraveling the complexities of nature and various technological applications.

A. RLC OSCILLATOR WITH LOCAL NEGATIVE IMPEDANCE
Earlier researchers have described that harmonic signals produced by oscillators can turn into a relaxation oscillation upon the change of a specific parameter.However, this type of transition is not observed in the well-known relaxation oscillator circuit depicted in Fig. 1, particularly when idealized circuit components are considered.To explore such behavior, we first focus on the simple RLC tank circuit shown in Fig. 3(a) and replace resistor R with a negative impedance converter (NIC), as shown in Fig. 3(b) in the second step.When energy is initially stored in the RLC tank at t = 0, the voltage V (t) will perform damped oscillations of the form when R is constant and R 2 > L/4C.The decay is determined by λ = 1/2RC.Amplitude A and phase φ depend on the initial conditions, whereas ω g represents the angular oscillation frequency.
Theoretically, sustained oscillations require G = 1/R = 0.However, even a slight deviation from this value causes either an exponential decay or a growing response.A growing response is associated with a negative value of R, and can only be generated by a dedicated active circuitry.In practical systems that maintain sustained oscillations, at least two energy-storage elements and one active nonlinear element featuring locally negative impedance (compensating for unavoidable losses of the components) are required.This locally negative impedance can be realized using a so-called negative impedance converter (NIC), shown in Fig. 3(b) which features an impedance of −R 1 between the shown terminals, provided that the output of the op-amp is not saturated.The op-amp output (assuming a rail-to-rail opamp) saturates at voltages . (For details on the NIC circuit, please see Appendix A and for further discussion see [12]).By attaching the NIC terminals of the circuit in Fig. 3(b) as a replacement for the resistor R in the tank circuit of Fig. 3(a) and applying Kirchhoff's current law at the node connecting the inductor, capacitor, and NIC, we obtain which becomes after differentiation with respect to time where G(V ) = dI /dV represents the differential conductance and I (V ) is the current-voltage characteristic of the NIC, as shown in Fig. 3(c).The value of G(V ) changes from −1/R 1 (domain A in Fig. 3(c)) to 1/R 2 (domains B) when saturation occurs.The differential equation (5) in dimensionless time where v is the dimensionless voltage The nonlinear system can be separated into two linear systems: System A for |v| < 1 and System B for |v| > 1, featuring two different eigenvalues resulting from a solution ansatz ∝ exp(λt): The associated individual solutions are given by with expansion coefficients K 1 to K 4 .As expected from the negative impedance, function v A (t) exponentially increases with time, whereas v B (t) decays.For a particular initial condition, for example, in terms of v and its time derivative at some time, depending on what range (e.g., A) is associated with these conditions, the respective expansion coefficients (e.g., K 1 and K 2 ) can be determined.The obtained solution is valid up to the time, where v crosses the border |v| = 1, where the FIGURE 5. Signal for nonlinear oscillator using ḠB = − ḠA = 5.The signal is composed by v A (red) and v B (blue) as given by ( 9) and (10).The dashed lines represent the single exponential functions in the solutions with expansion coefficients and the respective eigenvalues denoted.
respective other solution (B in our example) takes over.The coefficients (e.g., K 3 and K 4 ) can be determined from v and its derivative at that time.This solution is valid up to the next crossing of |v| = 1.Thus, the obtained partial solutions can be combined to obtain a complete solution for v(t). 1 Fig. 4 shows the voltage v over dimensionless time for ḠB = − ḠA = 0.3, 5, 10 and the small initial condition v(−50) = 10 −4 and v(−50) = 0.The oscillations are sinusoidal for small conductances (represented by eigenvalues that are approximately λ ≈ ±j), resulting in an angular oscillation frequency of ω ≈ 1/ √ LC (that is ≈ 1 using the scaled frequency ω = ω √ LC).However, for larger conductances, the oscillations become anharmonic, which can be explained by considering the eigenvalues in (7) and (8): For Ḡ2 A and Ḡ2 B much larger than 4, the eigenvalues are real-valued, very different, and approximately given by:  becomes increasingly anharmonic, which is a characteristic of relaxation oscillations.Both dimensionless conductances ḠA and ḠB become large as C tends to zero, leading to very sharp transitions and extended slow parts.As the relation between dimensionless and physical time is t = t/ √ LC, the slow part of the oscillation remains finite in physical time, even if C approaches zero.We refer to this situation, where the duration of the fast transition is negligible compared to the duration of the slow-part, as a fully developed relaxation oscillation.It is important to emphasize that if C is exactly zero, theoretically, only one energy-storing device remains in the circuit, and no oscillation can exist.However, in real circuits, some sort of spurious capacitance is always present, causing a fully developed relaxation oscillation.
Instead of pursuing deeper analyses of the associated differential equations (for more details see, e.g., [13]), an intuitive circuit analysis is performed to determine the characteristics of the fully developed relaxation oscillation with the resulting signals shown in Fig. 6, where it can be observed that the duration of fast transitions decreases to zero.
The circuit of the NIC in Fig. 3(b), loaded with an inductor, can also be viewed as a non-inverting amplifier with an amplification factor A 0 = 1 + R 2 /R 1 loaded with a series circuit consisting of L and R 2 .The node between L and R 2 is fed back to the positive input.From this viewpoint, the operating principle of this circuit can be derived conveniently.
Assuming that on turn-on, the output of the op-amp is saturated at the supply voltage V B , and there is no flux in the inductor (corresponding to I L,0 = 0), such that the potential at the non-inverting input is V B , which, for the moment, sustains the positively saturated output.Subsequently, as the magnetic field in the inductor builds up and the current starts to increase, the voltage V = V B − i L (t)R 2 across L decays exponentially until it falls below the level V = V B /A 0 required to saturate the op-amp output, which causes the output, and thus, V to decrease further.This positive feedback causes the output to quickly flip to the negative supply voltage −V B .Therefore, the voltage V jumps virtually instantly from V B /A 0 to −V B (2 − 1/A 0 ).After this rapid change, the current in the inductor continuously changes in the opposite direction with the same time constant as before.The associated waveforms are shown in Fig. 6.The circuit generates a square wave with a duty cycle of 50 % at the output of the op-amp.Consequently, the oscillation period is obtained using the charging equation of an LR series circuit where I L,0 and I L,∞ denote the initial and terminal inductor currents, respectively, analogous to (1).Based on the working principle discussed above, (13) can be transformed to give the period T , setting and I L (T /2) = −I L,0 , yielding Solving the nonlinear differential equation ( 6) for more complex nonlinearities is only possible using numerical methods for most problems.However, the insight gained from such an analysis is limited because important features of the oscillation (e.g., the period) cannot be directly obtained from the numerical model.Alternatively, intuitive circuit analysis reveals the operating principle and allows for the determination of the signals and oscillation periods.However, the result is accurate only for C ≈ 0 in the previously discussed example.The waveforms are not easily obtained when, for example, the amplifier features a soft saturation curve (as will be discussed later in this paper).The main aim of this work is to illustrate the use of a third method namely analysis of the fast-slow system representation of the equations, which is illustrated for the so-called van der Pol equation first and for various other circuits later.

B. THE FAST-SLOW SYSTEM REPRESENTATION AND THE VAN DER POL EQUATION
Equation ( 6) represents a special case of the so-called Liénard equation [14], which can be written as Although (6) allows for the exact solution discussed earlier, only a limited number of exact solutions to the general Liénard equation (15) are known (see, e.g., [15], [16], [17]).However, these known solutions are not directly related to circuit engineering problems; therefore, they are not discussed in this work.
The conditions for e(x) and f (x) in (15), for which unique and stable limit cycles exist, are obtained by considering the equivalent problem: with and applying the Liénard theorem [18], which states that under the assumptions • F and e are differentiable and odd in x, • xe(x) > 0 for x ̸ =0, • F(0) = 0 and F ′ (0) < 0, • F has a single positive zero at x = a, • F increases monotonically to infinity for x ≥ a as x → ∞ it follows that (15) has exactly one stable limit cycle.
Van der Pol derived a well-known equation based on studies of a triode vacuum tube oscillator circuit [1], as depicted in Fig. 2. Assuming a third-order polynomial transfer characteristic of the vacuum tube anode current, the dynamic equations can be transformed (see Appendix B) into the following dimensionless form [19] where we introduced a scaled time t (note that this was t in Appendix B), which will be re-scaled for convenience below.This equation represents an oscillating system with amplitude-dependent damping.For |x| < 1, the oscillation grows, whereas it decays for |x| > 1, yielding a stable oscillation after the settling period.However, because the damping is time-dependent, a pure harmonic oscillation can only be obtained for λ = 0.For larger values of λ, the time signal x exhibits fast transitions.The Liénard equation ( 15) is a generalization of the van der Pol equation (19) and can be represented by the so-called fast-slow system [20]: where two variables, i.e., a ''fast'' (x) and a ''slow'' (y) variable appear and replace the single variable in van der Pol's equation.If x and y are vectors of dimensions m and n the system is also termed (m, n)-fast-slow system; however, in this study, we consider only (1,1)-fast-slow systems.Parameter ε is positive and much smaller than 1.Characterizing the limit cycles of this generalized representation is difficult.
A reduced version of this problem, where g and h are polynomials, is known as Hilbert's 16th problem [21], which has only been partly solved to date.Even Poincaré abandoned the search for exact solutions in favor of studying qualitative features.Two important theorems are the uniqueness theorem, The state strictly follows the singular limit and approaches either one of the two points s 1 and s 2 , from which fast horizontal transitions occur also known as relaxations.
which states that trajectories in the phase diagram never intersect, and the Poincaré-Bendixson theorem [8], which states that if a trajectory is confined to a closed bounded region and there are no fixed points in the region, the trajectory must eventually approach a closed orbit.Throughout this paper, we use the convenient representation in ( 20) and ( 21) to describe various oscillator circuits by adapting the functions h and g to fit the respective problem.To obtain this fast-slow system representation for the van der Pol equation, the so-called Liénard transform (see [20]) is applied, where the new variable y is introduced as The scaled time in the Liénard equation is given by t = t/λ.Applying this to the van der Pol equation (19), we obtain: i.e., h = y−x 3 /3+x and g = −x.In the context of relaxation oscillations, the parameter ε = 1/λ 2 is sometimes referred to as the ''small parameter'', as the smallness of ε corresponds to fast signal transitions.To understand the reason for this, we first study the degenerate equation, where ε = 0: Setting ε = 0 corresponds to reducing one energy storage device in the circuit to zero.The cubic parabola in the first equation gives a strict relation between x and y in the phase diagram shown in Fig. 7 and is denoted here as the singular limit [20], whereas the second equation determines the dynamics in the phase diagram.The degenerate system does not represent a conventional differential equation in the sense that the initial conditions cannot be arbitrarily chosen [19].Instead, they must be points on the singular limit.For initial conditions, such as p 1 and p 2 (see Fig. 7), where x < 0, ẏ is positive, causing y to tend upward until point s 1 is reached.However, s 1 is not a stable point, because x is still negative at s 1 , therefore, y continues to tend upward (dy/dt > 0).Simultaneously, the state must remain on the curve defined by the first equation of the system.The only way to fulfill both conditions is by an instantaneous jump to the right branch of the curve, where x is positive.From here, point s 2 is approached (as now dy/dt < 0), where, similar as before, a jump to the left branch of the curve is induced.By eliminating y from ( 25) and ( 26), the time-derivative of x can be obtained as follows: Instantaneous jumps occur at points x(t) = ±1.These jumps appear nonphysical and are present only because of the oversimplification by setting ε to zero.When considering specific cases, it turns out that this oversimplification is commonly associated with neglecting a component or, more generally, a mechanism that requires the system to be of second order (otherwise, a point in the phase diagram cannot reverse its direction [8]).In the case of nonzero ε, the trajectories are allowed to depart from the singular limit.This deviation is equal to h in (20), and by moving ε to the right, the expression reveals that fast transients in x (i.e., fast horizontal movements in the phase diagram) can occur, particularly for small ε.
Returning to the consideration of the degenerate system ε = 0, points p 3 and p 4 are expected to tend toward s 1 or s 2 .However, the dashed part of the singular limit in Fig. 7 is unstable; therefore, for the initial conditions p 3 and p 4 , the state jumps to the stable parts of the curve.This can be understood by considering (27).When 0 < x < 1 then x moves in the positive x-direction and accelerates as it approaches 1.The stable parts are characterized by ∂h(x, y)/∂x < 0 and the unstable parts by ∂h(x, y)/∂x > 0. As pointed out in [19], these conditions can be explained kinematically in terms of the stability of the equilibrium positions of an auxiliary first-order equation by applying Liapunov's second method for stability [18].Relaxations occur at transition points s 1 and s 2 where ∂h(x, y)/∂x = 0.
(29) Fig. 8 shows examples for different ε values under various initial conditions.When ε is very small, fast horizontal transients occur, quickly bringing the state close to the singular limit.As the state approaches the singular limit up to a distance on the order of ε, it begins to follow the singular limit in close vicinity until it reaches the vicinity of s 1 in Fig. 7.At this point, a fast horizontal transition catapults the state against the positive arm of the singular limit, and subsequently follows the curve until it approaches point s 2 .A rapid transition to the negative leg of the singular limit closes the limit cycle.Fig. 8 shows that the fast horizontal transitions become ''smoother'' for larger ε, which is plausible when considering that the direction of the trajectory is dy/dx = −εx/h, i.e., the ratio of ( 26) and (28).
The oscillation frequency of the van der Pol oscillator can be estimated from the fast-slow system representation by assuming ε = 0 [7].align (21) characterizes the slow signal parts, which determine the oscillation period calculated by separation of variables: Here, dy is obtained from (25) and for the van der Pol oscillator this yields Therefore, the dimensionless period of oscillation is [8] However, the calculation becomes much more complicated for nonzero ε.Typically, asymptotic methods such as those presented in [22] must be employed.The result in (32) is considered to be a zeroth-order approximation of such an asymptotic expansion [19].The period of the van der Pol oscillator is known to be approximately 20% longer for ε = 0.01 [23], indicating that this approximation is coarse.However, in the following examples, the small parameter ε is introduced by considering spurious components, including wire inductances, the band limit of an op-amp, and the spurious capacitance of an LED.As a result, the actual ε is much smaller, and the period is accurately determined by the zeroth-order approximation.

II. VARIOUS RELAXATION OSCILLATORS
In the following sections, the transformation of the differential equations and the analysis of the fast-slow system representation, as shown in the example of the van der Pol oscillator, are applied to various oscillator circuits.These circuits include a saturating op-amp oscillator, the Schmitt trigger oscillator already introduced, and a popular simple voltage step-up converter known as the ''Joule thief''.All these oscillators share the characteristic that their principle of operation can be easily understood through inspection.The oscillation frequencies of these circuits can be determined by considering the charging and discharging of capacitors and inductors.However, although establishing a simple form of nonlinear ordinary differential equations (ODEs) for each circuit is most often feasible, merely inspecting these ODEs typically does not yield the signals or oscillation frequencies as obviously as intuitive circuit analysis does.Because advanced analyses rooted in nonlinear dynamics are required to obtain these insights, understanding circuits solely by their governing equations may appear unnecessarily complicated.
Nonetheless, the transformation to a fast-slow system representation, which is easily obtained in most cases, allows for insights beyond basic circuit analysis.This representation offers a powerful tool for exploring the behavior of oscillator circuits and gaining a deeper understanding of their dynamics and characteristics.

A. RELAXATION OSCILLATOR USING SATURATION 1) BASIC CIRCUIT ANALYSIS
In Fig. 9(a), an oscillator is shown, where the op-amp in the circuit is configured as a non-inverting amplifier with moderate amplification.The limited supply voltage causes the output voltage to saturate at a certain level.The transfer For engineers familiar with circuit theory, the principle of operation is easily revealed by, e.g., assuming that the capacitor is initially uncharged and the output initially provides a positive supply rail voltage (assuming that a rail-to-rail opamp is used for simplicity).As the capacitor charges, the voltage across the resistor decreases, as shown in Fig. 10.When the voltage across the resistor, which is also the input voltage of the op-amp, falls below the value required to keep the output in positive saturation, i.e., V B /A 0 , by virtue of the positive feedback of the then reducing output voltage, the output quickly swings to the negative rail −V B .Because the capacitor was close to fully charged the period T is derived as follows: The equality of the wave-forms in Fig. 6 and 10 is not coincidental.The step response of the resistor voltage of an RL circuit is identical to the voltage at the capacitor in an RC series connection for equal time constants.This circuit and the circuit in Fig. 3 may be considered unusual designs because op-amps with permissible input voltages outside the supply voltage range are rare.

2) FAST-SLOW SYSTEM REPRESENTATION
To obtain a fast-slow system representation, a small parameter ε must be introduced to yield a second-order differential equation.This can be achieved, for instance, by considering a small parasitic inductance L in series with the capacitor.The ODE for current I , with A(V ) denoting the transfer function 99458 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
of the amplifier (i.e., V O = A(V )) can be transformed into an ODE for the dimensionless input where A ′ = dA(v)/dv denotes the voltage-dependent differential small-signal amplification.A similar oscillator with an additional resistor R 2 in series to C and L is discussed in [24], where it is considered a negative resistance oscillator. 2 Indeed, (35) resembles that of an RLC circuit (when multiplied by R), featuring a nonlinear resistor R(1 For the nonsaturated op-amp, assuming A 0 > 1 and 4L ≪ CR 2 , the eigenvalues are approximately s 1 ≈ (A 0 −1)R/L and s 2 ≈ 0. Therefore, waveform v builds up exponentially and is dominated by v(t) = v(0) exp(s 1 t) for an initial condition v(0).Introducing a variable x = v and dimensionless time t = t/τ with τ = RC in (35) yields the following Liénard equation: The small parameter is ε = τ p /τ corresponds to the ratio of the time constants τ p = L/R to τ .Therefore, a small parameter ε also corresponds to a rapid voltage build-up during start-up.In Fig. 11 The exact solution can therefore be determined by piece-wise solving of two linear ODEs while maintaining continuity at the transition points, similar to the NIC example shown before.However, a graphical approach based on a fast-slow system representation will be pursued instead.By equating −x in (36) with dy/d t, the variable y is defined similar to the van der Pol equation: The integral in (37) is evaluated, yielding the associated fast-slow system representation as follows: For small spurious inductances L, ε is also small, resulting in relaxation oscillations with fast transitions.The equation is formulated in dimensionless time t = t/(RC), and the physical time axis scales with the time constant τ = RC. 2 It can be shown that the respective differential equation of the extended circuit can be transformed into (35) using the transformed quantities The singular limit, y = x − A(x) is shown in Fig. 12(b) for different amplification factors.Owing to op-amp saturation, a zig-zag curve results with clearly defined transition points at x t = ±V B /A 0 , as previously found by basic circuit analysis.
The oscillation period can also be determined directly for ε = 0 by using (30).For g(x, y) = −x = −(A(x) + y) (see (38)) , and considering the left slow-part in Fig. 12, where A(x) = −V B , the period in the dimensionless time frame is The integration bounds ±y p in (40) are determined by the local extrema of y at the singular limit.A(x) is either V B or −V B at the slow branches.The time constants obtained from this direct calculation and from the intuitive circuit analysis are equal when converted to the same time frame.For the case A 0 = 1 (40) yields a period of zero.However, it is apparent from (38) and (39) for A(x) = x, i.e., a harmonic oscillation with a period of T = 2π √ ε in dimensionless time and T = 2π √ LC in physical time is possible; however, it is not self-exciting.

3) SOFT SATURATION
The singular limit will look more similar to that of the van der Pol oscillator if a softer saturation function of the amplifier is considered.The oscillation period of such an oscillator depends sensitively on the saturation characteristics.Three different hypothetical sigmoid functions for the transfer function of the op-amp, including a hyperbolic tangent function, Gudermannian function, and simple algebraic function that saturate at ±1 featuring a slope of dA/dx = A 0 at x = 0, are considered: A gd (x) = 4 arctan (tanh (A 0 πx/4)) /π, (44) These functions are shown in Fig. 13(a) for various values of parameter A 0 .For hard saturation (sharp corners of the transfer function), the differential amplification dA(x)/dx jumps from A 0 to 0. However, for soft saturation, the transition points are not as obvious using basic circuit analysis, but they can be easily revealed by considering the fast-slow system representation (38) for ε = 0.According to (29), the transition points x t of the singular limit shown in Fig. 13(b) are given by the local extrema of y = x − A(x), associated with dy/dx = 0 (i.e., (38) for ε = 0) and hence Therefore, the transition points for soft saturation can be defined as, where the differential amplification drops to 1.
The limit cycle shown in Fig. 13(b) for the hyperbolic tangent function (43) resembles that of the van der Pol oscillator, and is excited by the same mechanism.Assuming the initial condition s 1 , the state approaches point s 2 from where a fast transition to s 4 is induced, followed by a slow approach to s 3 .The limit cycle is closed by another fast transition.From threshold −x t , an instantaneous jump to x p occurs.x p is therefore related to x t by y(−x t ) = y(x p ), i.e., −x t − A(−x t ) = x p − A(x p ).This can be rearranged for x p , yielding an attracting fixed-point iteration x p,n = A(x t )−x t +A(x p,n−1 ) as |A ′ (x)| < 1 on the slow branches.A suitable starting condition is, e.g., A(x p,0 ) = 1.To obtain the exact value for the oscillation period (assuming ε = 0), ( 30) can be used: The period in physical time is therefore given by The obtained period T is a function of A 0 and scales linearly with τ = RC.The numerical evaluation of the integral for various A 0 yields the results shown in Fig. 14.Surprisingly, the deviations of the periods increase with A 0 , although the saturation functions appear more similar as they approach the shape of hard saturation with increasing A 0 .The reason for this behavior lies in the differences of threshold voltages x t in (47) -(49).As the capacitor discharges, its voltage tends to zero, but relaxes as it reaches x t .Therefore, the period is particularly sensitive to x t for large A 0 .These results can be conveniently reproduced using LTspice simulation [25], as shown in Fig. 15.
The consideration of a spurious inductance in the circuit led to the identification of the small parameter ε.It can be shown that assuming low-pass behavior of the amplifier can be used equivalently to obtain a small ε parameter.This approach is used in the following example.

B. RELAXATION OSCILLATOR USING HYSTERESIS
In this example, the circuit discussed in the Introduction is analyzed using a fast-slow system representation.The typical interpretation of the circuit is that the capacitor voltage oscillates between the two thresholds ±V T of the Schmitt trigger.To introduce the small parameter ε in this example, a small time constant τ p is added to the output of the Schmitt 99460 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.trigger (see Fig. 16), which can be considered a first-order approximation of the slew limit or transition time of the opamp.The buffer is added for simplicity to decouple the RC circuit from the loads.The transfer function of the trigger is denoted as A(V ). Circuit analysis yields the following expressions:

1) FAST-SLOW SYSTEM REPRESENTATION
Introducing the dimensionless time t = t/τ where τ = RC and ε = τ p /τ with τ p = R p C p , and using the substitutions x = V O /V B and y = V /V B yields the fast-slow system representation The singular limit x = A(y) is shown in Fig. 17.For arbitrary initial conditions, the singular limit is quickly approached horizontally.Once in the vicinity of the singular limit, the state slowly follows the vertical segments with fast jumps between the transition points.
The formal calculation of the period from the dimensionless fast-slow system representation for ε = 0 yields the same result when converted to physical time, as was found by the circuit analysis in the Introduction, i.e.,

C. TRANSISTOR OSCILLATOR JOULE THIEF
The term Joule thief refers to a simple voltage step-up converter, which as far as we know, was first presented in [26].The circuit shown in Fig. 18, is designed to remain functional even at low battery voltages, enabling it to harvest (''steal'') the last remaining bits of energy from drained batteries that are no longer useful for conventional applications.One common application of the Joule-thief is to drive an LED flashlight.The circuit's appeal lies in its low part count and relatively simple transformer design, where a few turns on a small toroid are typically sufficient for operation.

1) BASIC CIRCUIT ANALYSIS
A basic circuit analysis is performed using idealized transistor and LED characteristics.The input (I C over V BE ) and output (I C over V CE ) characteristics of the transistor, as well as the current-voltage characteristics of the LED (I D over V D ), are approximated by abrupt functions, as shown in Fig. 19.
After connecting the battery to the circuit, a rapidly increasing small yet sufficiently large base current drives the transistor to saturation.Consequently, the voltage at the secondary side of the (tightly coupled) transformer jumps to the battery voltage V B minus the saturation voltage V sat of the transistor.Current I 2 increases linearly, as shown in Fig. 20, according to where L 2 is the inductance of the secondary winding.The primary current I 1 = −I B is small and its influence on the secondary side is neglected.The mutual inductance coupling the primary and secondary sides of the transformer is denoted by M , and the ratio of the transformer windings N 1,2 is denoted by a = N 1 /N 2 .The induced voltage V 1 = M İ2 = aV 2 on the primary side enhances the base current I B : driving the transistor further into saturation.However, when the constant base current, by virtue of the limited current gain β, is no longer sufficient to maintain the linearly increasing current I 2 , the transistor opens.At this point, the collector current reaches its maximum Î2 = βI B .The charging time T c is calculated using (58) by separation of variables and assuming I 2 (0) = 0, yielding: After the charging interval, the voltage V 2 inverts as inductance L 2 maintains the continuity of I .Simultaneously, the inverted primary voltage V 1 of the transformer off the transistor even more.During the discharging interval, the inductor drives the current through the LED with a forward voltage V LED assumed to be constant.Consequently, V 2 = V B − V LED is negative, leading to a linear decrease in I 2 according to L 2 dI 2 /dt = V 2 until the current almost vanishes.Hence, the discharging interval T d is given by The full period is T = T d + T c , and the duty cycle D can be calculated by The current delivered to the LED decreases linearly (due to the assumption of a constant V LED ) from Î2 to 0 during time T d (see Fig. 20) and is zero during the charging period T c .The RMS value of the current, denoted by ĪD , is calculated by This analysis demonstrates that the power delivered to the LED increases with higher values of a and β, and increases with a lower value of R. The most convenient method for adjusting LED power is therefore to vary R. The waveforms presented in Fig. 20 for the idealized circuit can also be simulated in LTspice by using the circuit shown in Fig. 21.The following transistor parameters are used for the simulation: β = 152, V sat = 0.28 V and V BE = 0.7 V.The LED's forward voltage is set to V LED = 2.1 V to achieve the best agreement with the actual circuit.Note that an abrupt change in the output characteristics of the transistor when leaving saturation results in a non-converging simulation.To address this issue and ensure numerical stability, a finite steepness parameter xi = 1µ is introduced, causing the collector current to increase from zero to the maximum value within a V CE = 1µV.In the simulation, a small capacitor C1 is incorporated into the idealized circuit to facilitate oscillations.However, its value can be as low as attofarads.
The results of this intuitive circuit analysis are approximate.In the actual circuit shown in Fig. 18, V sat , V LED , V BE , and β change over a cycle.This introduces additional complexities that are considered using the fast-slow system representation discussed below.

2) FAST-SLOW SYSTEM REPRESENTATION
The circuit in Fig. 18 can be analyzed in the phase diagram based on the fast-slow system representation where the variations of the transistor and LED parameters during one cycle can be considered, yielding a deeper insight into the operation of the circuit.To derive the fast-slow system representation, a small parasitic capacitance C is assumed to be in parallel to the LED.By circuit simulation, it was found that considering the LED voltage x = V D and the inductor current y = I 2 as the fast and slow variables yields a familiar limit cycle trajectory.Using Kirchhoff's voltage law at the primary and secondary sides of the transformer, and the current law at the collector and LED node, we can establish the following relations: Since the current I 1 is assumed negligible, the transformer equations reduce to The ratio a = M /L 2 equals the turns ratio of N 1 : N 2 of a tightly coupled transformer, and therefore Combining (65) with (68) yields which, together with (66), is a fast-slow system representation 3 for ε = C, x = V D and y = I 2 , when the collector current I C and the LED current I D are expressed in terms of V D and I 2 .The base-emitter voltage V BE depends on the base current I B (= −I 1 ) and the collector-emitter voltage V CE (= V D ).The LED current I D depends on the LED voltage V D , and as driving the transistor into saturation is a crucial part of the operation principle, the collector current I C must be considered as a function of the collector-emitter voltage and the base current, i.e.: Functions F, G, and T are nonlinear and can be obtained from circuit simulations.The output characteristics of the used transistor (BC547B) I C = T (V D , I B ) are depicted in Fig. 22 for constant base currents I B generated using the part model included in LTspice.Using (64) and (65) with V 1 = aV 2 and V BE = F(I B , V D ), an implicit relation between the base current I B and the LED voltage V D can be established as I B is therefore a function of the variable V D denoted by H , i.e., which is shown in Fig. 23 for various values of resistor R in the circuit depicted in Fig. 18.The solid lines in Fig.  are simulated, and the dashed lines are approximated using (59), assuming a constant V BE = 0.7 V.The bends upwards at low V D are associated with drops of V BE as the usually backward biased base-collector diode begins to conduct for V D < V BE .The characteristic voltage V K where the base current approaches zero, and the approximate base current for V D = 0 (i.e., I K ) are therefore given by Fig. 24 shows the collector current for various values of R under the constraint that the base current depends on V D (= x) through (74).The collector current I C can consequently be expressed as a function of the LED voltage alone, i.e., 66) and ( 69), the fast-slow system representation where x = V D , y = I 2 and ε = C is obtained.The singular limit of that system therefore is simply the sum of the collector (Fig. 24) and the LED current as a function of LED voltage x, i.e., y = T ′ (x) + G(x), giving a typical N-shaped curve shown in Fig. 25, which is needed for relaxation oscillations to occur.Here, T ′ (x) is shown for various values of R and G(x), i.e., the voltage-current relation of the LED, is shown for LEDs of three different colors following the PSice models in Appendix C. The forward voltages of the LED are all higher than V K (75) at which the collector current vanishes, such that transistor and LED conduct alternately.The critical point where relaxation to the right (i.e., a rapid transition) is induced is at the local maximum Î2 of I 2 , where the transistor leaves saturation at approximately x = V D ≈ 0.28 V (see Fig. 22).The base current is well-defined at this point (see Fig. 23), but the current gain β at the edge of saturation may depend sensitively on component variations.Consequently, the peak current driven through the LED, although deterministic in the simulation, may vary considerably in the experiment.Fig. 26 shows the waveforms for V D and I 2 and Fig. 27 for V BE and I B , respectively.Compared to the waveforms in Fig. 20 of the idealized circuit, it is observed that the forward voltage of the LED and the saturation voltage change considerably over time.The secondary current I 2 depletes to zero for V B = 1 V and a = 1.By increasing the battery voltage to, e.g., V B = 1.5 V, I 2 does not deplete, as can be observed from the singular limits shown in Fig. 28.Such a condition occurs when V K in (75) is larger than the forward voltage of the LED, which results when V B is increased or when a (the turns ratio) is lowered to a < 1.The minimum secondary current and the minimum voltage V D are nonzero in this case and are denoted by Ǐ2 and Vsat , respectively.The oscillation period is again calculated using (30): For the left slow branch, only the transistor (corresponding to T ′ ) conducts.An improved approximation to (60) for the charging time T c is obtained if the collector current I C over the collector-emitter voltage V D is linearly approximated between V D = 0 . . .Vsat , i.e., T ′ (x) = x Î2 / Vsat , which is justified when Fig. 24 is considered.The result is improved compared to that in (60) and reads: Here, Vsat denotes the peak saturation voltage, which can be obtained graphically from Fig. 25.The discharge time is calculated for the case when I 2 depletes using the diode current equation, considering the (typically considerable) series resistance R S (see Appendix C) where V th denotes the thermal voltage, and n the ideality factor.The discharge time T d follows from the integral which can be evaluated numerically when (81) is transformed for V D neglecting the 1 A good approximation for (eqn:sdfkjsdhf) can be obtained: To account for a non-depleting I 2 in case of V B = 1.5 V, the lower bound in the integral (83) can be replaced by Ǐ2 , which, however, does not consider that there is a decreasing collector current for V D < V K (see Fig. 25), which results in a dilation of the discharge time. 4n Tab. 1 the theoretical values and simulation results using the circuit in Fig. 29 are compared for a red LED and R = 2 kΩ at various supply voltages.The power efficiency, in terms of the ratio of the average power, delivered to the LED to the average power drawn from the battery (η in Tab. 1) assuming a lossless transformer, reaches its maximum of 91.5% for V B = 1.28 V.Under these optimum conditions, the collector current does not fully vanish with Ǐ2 = 7.7 mA.
In addition to V B and R, also the turns ratio a can be changed.Increasing a causes higher currents I K , and with it, higher LED peak currents, which also leads to lower required TABLE 1. Simulation results for a Joule thief circuit using a BC547B transistor with R = 2 kΩ and a red LED.The values in brackets were calculated using the respective equations given in the header of the table.minimum battery voltages, which are both desirable features of a Joule thief.The turns ratio a is therefore a parameter influencing the singular limit.Increasing the inductances while keeping the turns ratio constant influences therefore only the oscillation period, provided that the inductor core does not saturate.

III. SUMMARY AND CONCLUSION
The relaxation oscillator circuits in this work were analyzed using conventional circuit analysis and graphical analysis of the phase diagram based on the fast-slow system representation.In particular, it was shown that the graphical representation of the singular limit yieldes additional insights.For instance, one can estimate the oscillation onset for an arbitrary initial condition, determine more exact expressions for the oscillation periods, or find the parameter ranges for which relaxation oscillations are excited at all.Particularly, the singular limit reveals hysteresis behavior in oscillator circuits even if there is no obvious hysteresis present in a transfer function of the circuit's building blocks.
It should be emphasized that effects associated with relaxation oscillations can be much more faceted than those presented here.For instance, a modified van der Pol system, as discussed in [7], exhibits so-called ''Canard cycles'' for certain parameters, where fast horizontal transitions in the phase diagram are induced not at the local maxima of the singular limits.
The examples discussed in this paper are all of order two.Higher-dimensional systems, denoted as (m, n) fast-slow system representations may show chaotic behavior for certain parameters, where no closed orbits result [13], [27].
The mathematically more involved oscillators of fractional order [28] are often encountered when distributed systems, naturally represented by partial differential equations, are approximated using lumped elements yielding ordinary differential equations.An op-amp realization of an oscillator with fractional capacitance is discussed in, e.g., [29]).
The so-called two-stroke oscillations, known for some time [30], have recently been discussed for memristor circuits [31].The limit cycles of these circuits show one fast and one slow transition instead of two fast and two slow transitions.The examples which were discussed in this work are therefore also referred to as four-stroke oscillators.
In summary, even though there are myriad manifestations of relaxation oscillation phenomena in natural sciences and technology worth discussing, the authors believe that electrical engineers, who add the relatively simple technique presented in this paper to their toolbox, will benefit from the possibility of viewing oscillation phenomena from a different angle and tackling problems involving nonlinear differential equations more intuitively.

APPENDIX A THE NEGATIVE IMPEDANCE CONVERTER
In the following, the characteristics of the NIC circuit in Fig. 30 discussed in Sec.I are derived for the general case, where all resistors are potentially different.Assuming an op-amp featuring infinite open-loop amplification and whose output is rail-to-rail, the voltage V is also present across resistor R 1 by virtue of the negative feedback.The op-amp output voltage is therefore 99466 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.As the op-amp output voltage cannot exceed V B there is a threshold voltage for V at which the characteristics of the NIC change.In particular, if |V | > V B , the output of the op-amp will be pinned to ±V B .The input current of the NIC is therefore and the resulting differential resistances are Fig. 31 shows the corresponding I/V characteristics.In the case of R 2 = R 3 , as for the circuit in Fig. 3(b), it follows that the NIC features a negative resistance of −R 1 .Maximum current ±I max flows at the transition points at ±V T At this point also the output current from the op-amp is maximum

APPENDIX B VAN DER POL'S EQUATION
The following brief derivation of van der Pol's equation from the oscillator circuit shown in Fig. 2 follows coarsely the original paper [1].The anode current of the triode I a generally depends on the grid voltage V g and the anode voltage V a but can be approximately written as a function φ of a single variable, i.e., the ''lumped voltage'' given by V a + gV g , where g represents the ''voltage-ratio'' of the tube (also known as amplification factor).The grid current is neglected.
We only consider changes around the steady, though unstable, state (i.e., V a0 = V B , I a0 and V g0 = 0) and introduce the variables for the deviations v and i from these bias values: i = I a − I a0 = φ(V a0 + v + gV g ) − φ(V a0 ).( 96) The voltage at resistor V R is related to the battery voltage V B and anode voltage V a by As the (tightly coupled) transformer converts only AC signals with the voltage ratio −M /L, the grid voltage V g is zero-mean and therefore, V g = −vM /L and the associated variation in the anode current is with k = gM /L − 1. Differentiating the relation −i = i R + i C +i L , obtained using Kirchoff's law of currents, with respect to time and using equation (97), we obtain the differential equation Current i can be expanded into a Taylor series around V a0 (using the notation φ ′ = dφ(x)/dx etc.), i.e., To obtain the van der Pol equation, it is assumed that the transfer characteristic of the tridode can be approximated in the vicinity of the steady state by using the first-and third-order terms of (100) (see Fig. 32).Therefore, the time-derivative of i is approximately given by: Combining equations (99) and (101) and considering that in the steady state V a0 = V B (because V R is zero as there is no DC voltage drop across the inductor parallel to R) yields Introducing a dimensionless time t = t/ √ LC and multiplying by L yields By defining the dimensionless parameter λ as and introducing the transformed variable x finally, the van der Pol equation results: Self-excited oscillations occur for the conditions φ ′ (V a0 ) > 0, φ ′′′ (V a0 ) < 0 and λ > 0, i.e., gM /L > 1 + 1/(Rφ ′ (V a0 )).

APPENDIX C PSpice MODELS
The model for the bipolar transistor BC547B used in the simulation is that of the manufacturer NXP and can be found in the LTspice library provided by Analog Devices Inc.The LEDs in the simulation are conventional 5mm (T1 0.75) types with the data provided by the manufacturer OSRAM. red: .MODEL LO_541B-typ D IS=661.43E-24N=1.6455~RS=4.8592green: .MODEL LT_543C-typ D IS=10.000E-21N=1.6963~RS=3.6653white: .MODEL LW_541C-typ D IS=414.48E-15N=5 RS=22.499

FIGURE 1 .
FIGURE 1.The relaxation oscillator (a) uses an inverting op-amp Schmitt trigger and an RC circuit.The transfer characteristic of the Schmitt trigger is shown in (b) with the switching thresholds given byV T = V B R 1 /(R 1 + R 2 )(for a rail-to-rail op-amp).The capacitor voltage V oscillates between ±V T giving a square output signal V O .

FIGURE 3 .
FIGURE 3. a) A simple RLC parallel circuit.b) The negative impedance converter (NIC) replacing R in circuit (a) .c) The current over voltage relation of the NIC.

FIGURE 4 .
FIGURE 4.The NIC-loaded LC circuit performs approximately harmonic oscillations for small conductances, but the signal becomes anharmonic and shows the characteristics of relaxation oscillations.The time interval from −50 s to 0 s where the signals slowly build up, is not shown for a clearer representation of the developing steady-state signals.

Fig. 5
Fig.5shows the contributions v A and v B for ḠB = − ḠA = 5.The exponentially growing signal contribution v A is dominated by the large eigenvalue λ A+ = 4.79.As outlined above, the solution v B , composed of exponentially decaying contributions, takes over for |v| > 1. Somewhat unexpectedly, the signal increases further for the first moments due to the large negative eigenvalue λ B+ = −4.79before it slowly decays, dominated by the smaller eigenvalue λ B− = −0.21.Generally, that increasing dissimilarity of the magnitudes of the eigenvalues λ A− and λ A+ as well as that of λ B− and λ B+ , as Ḡ increases, causes a narrowing of the fast signal contribution v A while dilating the slow v B .The signal

FIGURE 6 .
FIGURE 6. Voltage signals for the circuit shown in Fig. 3 for V B = 15 V, L = 1 H, R 2 = 1 kΩ, C → 0 F, and A 0 = 5.To show the current I L on one common voltage axis it is multiplied by R 2 .

FIGURE 7 .
FIGURE 7.Phase diagrams for the fast-slow system representation of the van der Pol equation for different initial conditions p 1 to p 6 and ε = 0.The state strictly follows the singular limit and approaches either one of the two points s 1 and s 2 , from which fast horizontal transitions occur also known as relaxations.

FIGURE 8 .
FIGURE 8. Phase diagrams for the fast-slow system representation of the van der Pol equation for different ε and the initial conditions indicated by the colored dots.Fast relaxations occur for small ε and are associated with the horizontal segments in the phase diagram.

FIGURE 9 .
FIGURE 9. a) Relaxation oscillator using an amplifier with saturation.The small spurious inductance L is included to realize a small ε. b) The transfer function of the non-inverting amplifier is denoted by A and is shown for the two A 0 = 1 and 5.

FIGURE 10 .
FIGURE 10.Voltage waveforms for the circuit shown in Fig. 9(a) for V B = 15 V, C = 100 nF, R = 10 kΩ and A 0 = 5.The first peak of V is lower, as the capacitor is assumed empty at startup.
, the limit cycles for three different ε values are shown, along with their associated signals.The ideal saturation characteristic A(x) is shown in Fig. 12(a) for various A 0 values.A ′ (x) in (36) is a box function, taking the value A 0 in the linear range |x| < 1/A 0 and 0 when saturated.

FIGURE 11 .FIGURE 12 .
FIGURE 11.Phase diagrams for the oscillator using a saturating op-amp for supply voltage V B = 10 V and A 0 = 5 and different initial conditions and ε.The input voltage swings above the power supply which requires a careful selection of the op-amp.

FIGURE 13 .
FIGURE 13. a) Transfer functions (43) to (45) of smoothly saturating amplifiers and for various A 0 .b) The singular limits y = x − A ta (x) for the hyperbolic tangent and various A 0 .The transition points s 2 and s 4 and the signal maxima are s 1 and s 3 are indicated for A 0 = 3.The transition points lie on the dash-dotted line for all 1 ≤ A 0 , showing the high sensitivity to A 0 when it is close to 1.

FIGURE 15 .
FIGURE 15.LTspice schematic using the behavioral voltage source B1 to implement the tanh saturation characteristic in (43) for A 0 = 100.Replace the expression V=tanh(V(x) *A0) by one of the functions in the comments (blue text) for the other characteristics.The period of oscillation T is found in the log-file.The current source I1 provides a short current pulse at t = 0 to initiate the oscillation.

FIGURE 16 .
FIGURE16.Relaxation oscillator using a Schmitt trigger from Fig.1where the small spurious time-constant τ = R P C p is included to realize a small ε.The transfer function of the amplifier is denoted by A.

FIGURE 17 .
FIGURE 17. Phase diagrams for the Schmitt trigger oscillator for dimensionless threshold voltage V T /V B = 0.66 and different initial conditions and ε.The spurious output R p C p causes the expected slower settling of the dimensionless output voltage x.

FIGURE 19 .
FIGURE 19.Idealized input (a) and output (b) characteristics of transistor and LED (c).I C = βI B for V CE > V sat .

FIGURE 21 .
FIGURE 21.An idealized circuit simulation using LTspice utilizing a behavioral current source (xi introduces a high, but finite steepness of the saturation characteristic required for numerical stability).

FIGURE 22 .
FIGURE 22.The output characteristic I C over V CE (= V D ) for a BC547B transistor for constant base currents I B , simulated using the part model provided by Analog Devices, see Appendix C.

FIGURE 23 .FIGURE 24 .
FIGURE 23.The simulated relation between base current and collector-emitter voltage I B = H(V D ) following (73) for a supply voltage V B = 1 V and a 1:1 transformer, i.e. a = 1.

FIGURE 25 .
FIGURE 25.Singular limit for various values of R and three different LED types and V B = 1 V and a = 1.A limit cycle is shown for a red LED and a base resistor R = 2 kΩ.

FIGURE 28 .
FIGURE 28.Limit cycles of the Joule thief using a BC547B with red LED and R = 2 kΩ for four supply voltages V B and a = 1.In cases where V K in (75) is larger than the forward voltage of the LED, the current does not fully deplete.

FIGURE 29 .
FIGURE 29.The LTspice circuit used to determine the values in Tab. 1.The results of the.meas directives are listed in the log-file of the simulation.

FIGURE 31 .
FIGURE 31.Current I over voltage V at the input terminal of the NIC circuit (black) and op-amp output voltage (blue).

FIGURE 32 .
FIGURE 32.Assumed characteristic of the anode current I a in equation (94) of the triode over the anode voltage V a when the grid voltage V g is expressed in terms of the anode voltage.