A nonlinear RLC resonator is investigated experimentally and numerically using bifurcation analysis. The nonlinearity is due to the parallel combination of a semiconductor rectifier diode and a fixed capacitor. The diode’s junction capacitance, diffusion capacitance, and DC current–voltage relation each contribute to the nonlinearity. The closely related RL-diode resonator has been of interest for many years since its demonstration of period-doubling cascades to chaos. In this study, a direct comparison is made of dynamical regime maps produced from simulations and circuit measurements. The maps show the variety of limit cycles, their bifurcations, and regions of chaos over the 2D parameter space of the source voltage’s frequency and amplitude. The similar structures of the simulated and experimental maps suggest that the diode models commonly used in circuit simulators (e.g., SPICE) work well in bifurcation analyses, successfully predicting complex and chaotic dynamics detected in the circuit. These results may be useful for applications of varactor-loaded split ring resonators.

Behavior of the resistor–inductor–diode oscillator is a popular example of complexity and chaos. The nonlinearity of the semi-conductor junction is accounted for by well-known models used regularly in circuit simulators, such as SPICE.1 This nonlinear oscillator is also the basis for nonlinear split ring resonators, which are of interest due to interesting magnetic properties produced when used to form arrays in metamaterials. Here, bifurcation analysis and Lyapunov exponents are used to calculate a dynamical regime map in the 2D parameter space of the source frequency and amplitude, which is compared to a regime map created from circuit measurements. Good agreement is found on the structure of the maps as indicated by the types of limit cycles, their bifurcations, and period-doubling cascades to chaos. The results suggest that these types of bifurcation analyses may be useful for developing applications of the nonlinear split ring resonators in metamaterials.

The series RLC circuit is one of the most common examples of the damped harmonic oscillator. In the numerical and experimental study presented here, nonlinearity is introduced into the resonator by including a diode in parallel with the capacitance. Semiconductor diodes have nonlinear current–voltage relationships, which are accounted for by mathematical models used in the well-known SPICE modeling tool for circuit design.1 The nonlinearities cause the circuit oscillations to have complex responses to a sinusoidal driving voltage. In this work, bifurcation analysis and Lyapunov exponents are used to find the limit cycles and chaotic attractors of this RLC-diode circuit and their dependence on the frequency and amplitude of the driving source. These simulation results are then compared to circuit measurements. The goal here is to determine how well the bifurcation and Lyapunov exponent analysis can predict actual circuit behavior over a wide range of parameter space.

The circuit considered here is essentially the same as the resistor–inductor–diode (RLD) circuit used earlier to demonstrate experimentally the period-doubling cascades to chaos, which can be induced by varying parameters of the sinusoidal voltage source.2,3 Those experiments also demonstrated the spacing convergence of the period-doubling bifurcation parameter values predicted by Feigenbaum’s theory on the universal behavior of period-doubling systems.4,5 More studies soon followed. Computer calculations based on the RLD circuit analysis predicted the period-doubling cascade to chaos.6 Intermittency routes to chaos were demonstrated.7 Comparison of model calculations and experimental measurements was made using bifurcation diagrams, return maps, and Poincaré plots.8,9 A simplified diode model was used to calculate a 2-dim bifurcation regime map, which agreed well with experimentally observed bifurcations.10 

The good agreement of prediction and experiment demonstrated in the studies mentioned above resulted in the RLD circuit becoming a popular example of a chaotic system. A good example relevant to this paper is the comparison of numerical predictions and circuit measurements by Carroll and Pecora showing good agreement on the location of the first period-doubling in 2-dim parameter space.11 

In the years since the above studies, bifurcation analysis tools have improved and computing power has increased. It is somewhat surprising that more use has not been made of these tools to make detailed predictions of the dynamical regimes of the RLD circuit by using the diode’s well-established mathematical models. To the best knowledge of the author, a direct comparison of numerically calculated and experimentally measured bifurcation regime maps, delineating the various limit cycles over 2-dim parameter space, has not been done for the RLD circuit.

In the past couple of decades, there has been increasing interest in exploiting nonlinear RLC circuits in metasurfaces for their interesting magnetic properties.12 In these applications, the RLC circuit is a split ring resonator (SRR) with nonlinearity introduced into the capacitance, often as a semiconductor junction,13–16 thereby making the SRR a RLD. The SRRs are typically small (dimensions of 10s of μm) and are driven by the electromotive force induced by electromagnetic waves in the GHz or THz range. The system considered here is conceptually and mathematically the same as a SRR. Only the size and frequency differ. Here, a function generator drives a circuit constructed from common electronic components: a resistor, an inductor, and the parallel combination of a fixed capacitor and a 1N4003 rectifier diode. The diode’s semiconductor junction is the source of the nonlinearity, just as it is for those SRRs, which incorporate varactor diodes.

The early studies of the RLD were primarily motivated by scholarly interest in understanding the nature of their complex and chaotic behavior. The more recent studies of arrays of the nonlinear SRRs used in metasurfaces are aimed at applications and, therefore, may benefit from improved accuracy of mathematical models of the nonlinear SRR.

The goal of this work is to compare simulations and circuit measurements of the RLC-diode over the 2-dim parameter space of the source’s amplitude and frequency. The simulations use bifurcation tools and Lyapunov exponent calculations to discover the limit cycles, their stability, and chaotic regions over the 2-dim parameter space. The model used for the diode accounts for junction capacitance, diffusion capacitance, and the nonlinear current–voltage relation by using well-established models. Circuit measurements detecting period-doubling, hysteresis, and chaotic-like behavior are made using a standard digital function generator and an oscilloscope. The work here refers to the RLC-diode circuit because the diode is in parallel with a fixed capacitor, whereas the RLD circuit generally does not include a fixed capacitor. The fixed capacitor allows for additional flexibility in lowering the resonance frequency of the circuit but does not change the essential nature of comparing the predictions and measurements of complex behaviors caused by the presence of the diode.

The presentation style here attempts to be accessible to those who are interested in the RLD circuit but are not particularly familiar with bifurcation analysis methods and terminology. Therefore, descriptions of the bifurcation methods are restricted to relevance to this work and may lack more general mathematical completeness and rigor.

The system of interest shown in Fig. 1 is a sinusoidally driven series RLC-diode circuit with a voltage-dependent capacitance due to the diode. We investigate how the dynamics of the circuit oscillations depend on the amplitude and frequency of the driving source. The focus is on using bifurcation analysis methods to identify the limit cycles and determine their stability.

FIG. 1.

Driven nonlinear RLC-diode circuit where a diode provides the nonlinearity. The semiconductor diode is modeled by the parallel combination of junction capacitance (Cj), diffusion capacitance (Cdf), and DC diode (the diode symbol). A fixed capacitance in parallel with the diode is also included. R represents the sum of explicit resistance and the inductor’s intrinsic resistance.

FIG. 1.

Driven nonlinear RLC-diode circuit where a diode provides the nonlinearity. The semiconductor diode is modeled by the parallel combination of junction capacitance (Cj), diffusion capacitance (Cdf), and DC diode (the diode symbol). A fixed capacitance in parallel with the diode is also included. R represents the sum of explicit resistance and the inductor’s intrinsic resistance.

Close modal

The capacitance of the RLC-diode resonator is the result of the parallel combination of the diode and the fixed capacitance (Cfx in the figure), which includes stray capacitance. The mathematical model used for the diode is the parallel combination of junction capacitance Cj, diffusion capacitance Cdf, and the diode’s standard DC current–voltage relation. All three of these terms are voltage dependent and, therefore, contribute to the nonlinearity of the system. The well-known models of semiconductor diodes used in the SPICE circuit simulator are used to account for these three effects.

For a voltage source V s = A cos ( Ω t ) driving the RLC-diode resonator, summing voltages around the circuit give
V s = L d I d t + I R + V d ,
(1)
where V d is the voltage across the parallel combination of the diode and fixed capacitor.
Current I from the voltage source splits, supplying current d Q / d t to the junction capacitance, current I d f to the diffusion capacitance, current I f x to the fixed capacitor, and current I d to the DC diode, giving the equation
I = d Q d t + d Q d f d t + d Q f x d t + I d .
(2)
Source current I and junction capacitance charge Q are the two variables chosen to describe the system. Therefore, voltages and other charges must be expressed in terms of I and Q.
I d f and I f x are expressed in terms of Q by using
d Q i d t = d Q i d V d d V d d t = C i d V d d t = C i d V d d Q d Q d t ,
where incremental capacitance C i = d Q i / d V d is used and subscript i = ( d f , f x ) for the diffusion and fixed capacitances. Solving Eq. (2) for the current into the junction capacitance gives
d Q d t = I I d 1 + ( C d f ( V d ) + C f x ) d V d d Q .
(3)
The  Appendix shows that the diode model allows V d and I d to be expressed in terms of Q. Thus, Eqs. (1) and (3) are the foundation for modeling the circuit in Fig. 1 in terms of the source current I and the diode junction capacitance charge Q.
Equations (1) and (3) are converted to a dimensionless form by defining dimensionless variables,
q Q Q 0 , i I V J L C J 0 + C f x , v V V J , t Ω 0 t ,
(4)
where V J and C J 0 are junction capacitance parameters (see the  Appendix). Q 0 = V J C J 0 is a naturally defined junction charge parameter and Ω 0 = 1 / L ( C J 0 + C f x ) is the nominal resonant frequency. Renaming t as t results in
d q d t = ( i i d ( q ) ) C J 0 + C f x C J 0 + ( C d f ( q ) + C f x ) d v d ( q ) d q ,
(5a)
d i d t = v d ( q ) σ i + μ cos ( ω t ) ,
(5b)
where t is now dimensionless. The dependence of v d on the variables ( q , i ) is required in Eq. (5). The  Appendix shows that
v d ( q ) = { 1 ( 1 ( 1 M ) q ) 1 1 M for  q < q t , F C + ( 1 F C ) M ( q q t ) for  q > q t
(6)
and
d v c a p ( q ) d q = { ( 1 ( 1 M ) q ) M 1 M for  q < q t , ( 1 F C ) M for  q > q t ,
(7)
where
q t = 1 ( 1 F C ) ( 1 M ) 1 M .
(8)
M and F C are two more junction capacitance parameters. The case for M = 1 is in the  Appendix.
The diode’s DC current–voltage relation requires two parameters, reverse saturation current I S, and emission coefficient N. The dimensionless relation is
i d ( q ) = i S ( e 40 V j v d ( q ) / N 1 ) ,
(9)
where Eq. (4) is used to obtain the dimensionless i S from I S. The diffusion capacitance is
C d f ( q ) = 40 T T × I S N e 40 V J v d ( q ) / N ,
(10)
where T T is the parameter for forward transit time.

We use parameter values appropriate for the 1N4001-1N4004 family of rectifier diodes: I S = 70 pA, N = 1.4, C J 0 = 33 pf, V J = 0.35 V, M = 0.45, F C = 0.5, and T T = 5 μs.

Remaining parameters in Eq. (5) are
σ = R C J 0 + C f x L , μ = A V J , ω = Ω Ω 0 .
(11)
In this work, ω and μ are the parameters of interest to vary in the bifurcation analysis. Equation (5) is the non-autonomous version of the equations for our driven RLC-diode circuit.
Some of the bifurcation analysis methods are more easily suited to autonomous systems. Therefore, Eq. (5) is converted to an autonomous form by including two additional variables u ( t ) and v ( t ). The differential equations for these variables are constructed so that u and v are sinusoidal with frequency ω and unity amplitude.17,18 The autonomous version of system (5) is then
d q d t = ( i i d ( q ) ) C J 0 + C f x C J 0 + ( C d f ( q ) + C f x ) d v d ( q ) d q ,
(12a)
d i d t = v d ( q ) σ i + μ u ,
(12b)
d u d t = u ( 1 u 2 v 2 ) ω v ,
(12c)
d v d t = v ( 1 u 2 v 2 ) + ω u .
(12d)
Initial conditions should satisfy u 2 + v 2 = 1. The terms containing 1 u 2 v 2 ensure that the amplitudes of u and v remain at 1. The sinusoidal behavior of u and v allows for periodic behavior of all four variables. This addition of two variables is in contrast to converting to an autonomous system by adding a single variable u ( t ) = t, which results in one fewer equation, but has the disadvantage of not allowing periodic behavior of all the variables.

The existence of the sinusoidal driving term ensures that there will be no fixed points (except for the trivial case of μ or ω = 0) and that the main limit cycle is period-1, meaning an oscillation with the same period as the sinusoidal source.

Bifurcation analysis of Eq. (12) was done using XPPAUT17 and AUTO-07p.18 These software tools are used to find the limit cycles (LC), their bifurcation points, and 2-dim parameter space regime maps showing the types of LCs. Time-series integration was done using an adaptive stepsize fourth-order Runge–Kutta routine. Lyapunov exponents were calculated using the method of Wolf et al.19 The bifurcation analysis consists of generating bifurcation continuation plots and the 2-dim regime maps. The bifurcations for Eq. (12) are period-doublings (PD) and limit points (LPs). Here, we describe the bifurcation points and their relation to the bifurcation continuation plots. It is worth noting that while the diode’s mathematical model used here is the same as the model in the SPICE simulator, the SPICE software itself was not used for any of the numerical simulations.

PD points occur on a bifurcation continuation plot at parameter values where a stable LC becomes unstable and a LC with double the period is created. The main LC of Eq. (12) is period-1, having the same period as the source. This LC exists as either stable or unstable over the entire parameter space. Figure 2 shows the ω-continuation bifurcation plot for the main LC. Six PD points are seen, where the LC changes stability. The solid (dashed) line indicates stable (unstable) LC. The plot shows how i-max (the maximum value of variable i’s oscillation) depends on parameter ω. LPs (also called saddle-node bifurcations) occur on a bifurcation continuation plot at parameter values where a stable LC branch and an unstable LC branch extending together in parameter space reach an extreme parameter value (the LP) where they collide and disappear. Another description is that a LC branch reverses its direction along the parameter-axis at the LP, with a concomitant change in stability.

FIG. 2.

Numerical ω-continuation bifurcation plot showing the maximum value of i for the LC. μ = 2 and σ = 0.013. The solid (dashed) line indicates stable (unstable) LC. PD and LP bifurcation points are indicated. Insets show coexisting LC time-series of q (orange) and i (blue) at ω = 0.372.

FIG. 2.

Numerical ω-continuation bifurcation plot showing the maximum value of i for the LC. μ = 2 and σ = 0.013. The solid (dashed) line indicates stable (unstable) LC. PD and LP bifurcation points are indicated. Insets show coexisting LC time-series of q (orange) and i (blue) at ω = 0.372.

Close modal

The two LPs in Fig. 2 create a region of hysteresis between their locations at ω = 0.316 and 0.376, where there is coexistence of large and small amplitude LCs as demonstrated by the time-series insets. The large amplitude LC is stable over the hysteresis region, while the small amplitude LC becomes unstable for ω below 0.369 where it undergoes period-doubling.

The LP and PD bifurcation points detected in the parameter continuation plots are the starting points for tracing their paths in the 2-dim ω μ parameter space. The paths are the borders for different dynamical regimes. Creation of these dynamical regime maps is a goal of this study because of their usefulness in presenting how the system Eq. (12) behaves throughout the ( ω , μ ) parameter space.

ω-continuation bifurcation plots of stable period-doubled LCs are used to find higher order period-doubling points. Those points are then traced in the 2-dim ω μ parameter space to find their borders allowing determination of period-doubling cascade regions.

The Lyapunov exponents (LEs) are useful for detecting regions with no stable LCs, where the stable oscillatory states are either quasiperiodic or chaotic. Borders for chaotic regions are not detected by the bifurcations tools but are easily identified in LE plots as parameter values where the maximal LE changes from zero (indicating stable LC) to positive (indicating chaos).

Circuits were constructed based on Fig. 1. The inductor was L = 10 mH with 23 Ω of intrinsic resistance. The diode was 1N4003, a standard rectifier. Fifteen of these diodes were used, and in-house measurements were made of the zero-bias junction capacitance, the reverse breakdown voltage, and the reverse-recovery time. Junction capacitance ranged from 23 to 85 pf, breakdown voltage was 300–1400 V, and reverse recovery time was 8–32  μs. All 15 diodes were tested in RLC-diode circuits to check the consistency of the dynamical regime results when using different diodes and to determine the sensitivity to the diode properties.

The capacitor in parallel with the diode was 68 pf. The fixed capacitance C f x in Fig. 1 accounts for both the 68 pf capacitor and stray capacitance. Stray capacitance was determined by measuring the resonant frequency for a small signal diode with a nominal 1 pf junction capacitance, yielding stray capacitance of 24 pf. The resistor was 100 Ω making the total series resistance R = 123 Ω. Using 123 Ω with C f x = 92 pf and C J 0 = 33 pf (for the particular diode used to create the experimental results regime map) in Eq. (11) gives σ = 0.013 for use in simulations.

The resonance frequency was determined using a small amplitude (25 mV) source. Figure 3 shows the frequency response of the diode voltage amplitude for source amplitudes of 25 and 100 mV. The voltage amplitude has been normalized by the source amplitude. The 25 mV source plot shows a typical resonance response, with a resonant frequency of 140 kHz. The maximum forward bias is about 12 × 25 mV = 0.3 V. The 100 mV source plot shows the distorted resonance response, which is due to the significantly higher diode current at the maximum forward bias of 5 × 100 mV = 0.5 V.

FIG. 3.

RLC-diode circuit’s frequency response. Normalized diode amplitude for 25 and 100 mV source amplitudes. R = 123 Ω, L = 10 mH, and 1N4003 diode in parallel with a 68 pf capacitor and 24 pf stray capacitance.

FIG. 3.

RLC-diode circuit’s frequency response. Normalized diode amplitude for 25 and 100 mV source amplitudes. R = 123 Ω, L = 10 mH, and 1N4003 diode in parallel with a 68 pf capacitor and 24 pf stray capacitance.

Close modal

Figure 4 is the numerical dynamical regime map showing the limit cycles and chaotic oscillations of the RLC-diode circuit in the 2-dim ω μ parameter space. The map consists of LP and PD lines. Traversing the map at μ = 2, there are six crossings of first PD lines and two crossings of LP lines, corresponding to the six PD points and two LP points in the ω-continuation of Fig. 2.

FIG. 4.

Numerical ω μ regime map for σ = 0.013, C f x = 92 pf. P1 (period-1), P2, and P4 LC regions are indicated. Hatching indicates higher order doublings and chaos. Frequency ω and source amplitude μ are dimensionless.

FIG. 4.

Numerical ω μ regime map for σ = 0.013, C f x = 92 pf. P1 (period-1), P2, and P4 LC regions are indicated. Hatching indicates higher order doublings and chaos. Frequency ω and source amplitude μ are dimensionless.

Close modal

A region of hysteresis of the period-1 (P1) main LC is indicated by the LP lines. Inside this narrow region, two versions of the main LC coexist, which causes a jump in the LC amplitude during a sweep of the frequency when a LP is encountered causing a sudden switch from one version to the other. Figure 2 shows the hysteresis and the time-series of the two coexisting P1 LCs for μ = 2.

There are two regions of period doubling of the P1 LC, a narrow region at a low frequency and a much larger region, which occupies a substantial amount of the map. The larger PD region contains two regions of higher order period-doubling cascades. The second and third doubling lines are shown. Period-doubling cascades to chaos are common when the spacing between higher order doublings decreases, as is the case here. The third period-doubling line can often be used as an approximation for the doubling cascade’s onset of chaos.

Figure 5 shows a bifurcation diagram corresponding to a traverse across the Fig. 4 regime map at μ = 4. The period-doubling cascades to regions of presumed chaotic behavior agree with those seen in the regime map. Confirmation of chaos by calculation of Lyapunov exponents is shown below.

FIG. 5.

Bifurcation diagram for increasing source frequency ω with σ = 0.013 , μ = 4. First, second, and third period-doubling points are indicated.

FIG. 5.

Bifurcation diagram for increasing source frequency ω with σ = 0.013 , μ = 4. First, second, and third period-doubling points are indicated.

Close modal

LE calculations detect parameter ranges with chaotic behavior and are, therefore, complementary to the parameter regime maps created by bifurcation analysis since those methods do not detect chaos. In LE plots for an autonomous system, stable LC is indicated by the maximal LE being zero, whereas chaotic behavior has a positive maximal LE. The LE plot for source amplitude μ = 4 in Fig. 6 shows that there are two regions where the maximal LE is positive (where the blue LE exceeds the red LE), near ω = 0.5, and a wider region around ω = 1.35. These two regions confirm the period-doubling cascades to chaos suggested by the nested higher order doublings in Fig. 4. Many such occurrences are labeled in Fig. 6 based on their agreement with where PD lines are crossed at μ = 4 in Fig. 4.

FIG. 6.

Three maximal LEs for increasing source frequency ω with σ = 0.013 , μ = 4. Chaos occurs where the blue LE is positive. First, second, and third period-doubling points are indicated.

FIG. 6.

Three maximal LEs for increasing source frequency ω with σ = 0.013 , μ = 4. Chaos occurs where the blue LE is positive. First, second, and third period-doubling points are indicated.

Close modal

Information about PD points can also be indicated on the LE plot. When a changing parameter causes a stable LC to undergo period doubling, the maximal LE remains at zero throughout the process. However, the second maximal LE increases until it reaches zero at the PD bifurcation point, and then it returns to increasingly negative numbers.

Figure 7 shows the simulated phase plot of the source voltage vs the diode voltage at ω = 1.25 , μ = 3.95 inside the right-side period-doubling cascade to chaos in Fig. 4. The phase plot shows the “filling in” of regions of the phase space characteristic of chaos.

FIG. 7.

Numerical phase plot for μ = 3.95 , ω = 1.25 , σ = 0.013, located inside the period-doubling cascade.

FIG. 7.

Numerical phase plot for μ = 3.95 , ω = 1.25 , σ = 0.013, located inside the period-doubling cascade.

Close modal

The LE plot gives information only for the dynamical state whose basin of attraction has captured the system. This is in contrast to the regime maps, which can show bifurcation lines for coexisting states, regardless of their basins of attraction. This means that bifurcation lines on the regime maps might not be indicated in a LE plot. An example is the right-side second PD of the left doubling cascade in Fig. 4, which shows this PD at ω = 0.66 for μ = 4. The third PD is very close to the second, suggesting that the cascade to chaos should also occur close to ω = 0.66. However, Fig. 6 shows that chaos extends to only 0.61. The existence of the two third PD lines that overlap indicates a complicated coexistence of higher order doublings and two separate cascades to chaos. Presumably, the basin of attraction of a stable period-doubled LC is able to capture the chaotic state at ω = 0.61.

PD points were detected as a function of source frequency and amplitude by varying a source frequency at each source amplitude while monitoring the diode voltage waveform. The change in the waveform indicating period doubling was easily recognized. Figure 8 shows the results for the circuit in Fig. 1 using R = 123 Ω ( 100 Ω resistor added to inductor’s 23 Ω) and 68 pf capacitor (added to 24 pf stray capacitance giving C f i x = 92 pf), giving resonance frequency 140 kHz. The 10 mH inductor is known to have nonideal properties (parasitic capacitance), which become problematic for frequencies above 300 kHz.

FIG. 8.

Experimental ω-“source amplitude” regime map showing PD borders and hysteresis induced amplitude jumps. Measurements from 1N4003 diode, L = 10 mH, and R = 123 Ω. f 0 = 140 kHz.

FIG. 8.

Experimental ω-“source amplitude” regime map showing PD borders and hysteresis induced amplitude jumps. Measurements from 1N4003 diode, L = 10 mH, and R = 123 Ω. f 0 = 140 kHz.

Close modal

Hysteresis was detected by noting the frequencies at which there was a large jump in the amplitude of the diode voltage LC. The location of the hysteresis for the circuit agrees well with the region in the numerical simulation, although it is much narrower for the circuit measurements.

The overall structure of the regime maps in Figs. 4 and 8 for simulations and circuit data has much in common. It is also informative to compare the time series of the simulations and the circuit measurements. Figure 9 shows good agreement of the LC period doublings, which occur when a source amplitude is increased at the nominal resonance frequency ω = 1 in both the numerical simulations of Fig. 4 and the circuit measurements in Fig. 8. The bottom row shows that for the small source amplitude (blue), the diode voltage LC (orange) has the same period as the source and is nearly sinusoidal. The middle row shows the time series for the period-doubled P2 LC. The general shape of the doubled LC is very similar for simulation and circuit measurement, and the LC’s period clearly requires two oscillations of the source. The top row shows the second doubling to P4 LC, which means that the LC requires four oscillations of the source. Again, the general shapes of simulation and measurement have good agreement.

FIG. 9.

Time-series simulations (left column) and circuit measurements (right column) for 1N4003, L = 10 mH, and R = 123 Ω showing period-doublings as the source amplitude increases for f 0 = 140 kHz. Blue (orange) is the source (diode) voltage. The bottom row is main LC, period-1. The middle row is period-2 LC after first PD. The top row is period-4 LC after second PD.

FIG. 9.

Time-series simulations (left column) and circuit measurements (right column) for 1N4003, L = 10 mH, and R = 123 Ω showing period-doublings as the source amplitude increases for f 0 = 140 kHz. Blue (orange) is the source (diode) voltage. The bottom row is main LC, period-1. The middle row is period-2 LC after first PD. The top row is period-4 LC after second PD.

Close modal

The one large discrepancy between simulations and circuit measurements is the behavior at high frequency and a low source amplitude where simulations show the P1 period-doubling region extending sharply down at twice the resonant frequency.

The existence of chaotic behavior in the circuit was investigated by continuing to increase the source amplitude beyond that used in Fig. 9 so as to complete the period-doubling cascade to chaos presumed to exist inside the nested doublings. At ω = 1 , ( f = 140 kHz ), Fig. 8 suggests that a source amplitude just above 2 V is likely beyond the onset of chaos. The oscilloscope phase plots in Fig. 10 show the behavior with a 1.85 V source, between the second and third doublings, prior to chaos (left panel), and with a 2.2 V source (right panel). The left panel shows the expected period-4 LC created by two period doublings, and the right panel shows behavior whose appearance is consistent with chaos. There is strong similarity between the right panel and the simulated chaotic phase plot in Fig. 7.

FIG. 10.

Oscilloscope screenshots of phase plots of measured diode voltage vs source voltage. (a) is period-4 LC with a 1.85 V source amplitude, and (b) is likely chaos with a 2.2 V source amplitude. f = 140 kHz.

FIG. 10.

Oscilloscope screenshots of phase plots of measured diode voltage vs source voltage. (a) is period-4 LC with a 1.85 V source amplitude, and (b) is likely chaos with a 2.2 V source amplitude. f = 140 kHz.

Close modal

As mentioned in Sec. IV, 15 1N4003 diodes were examined, which had a wide range of measured diode properties. All the diodes reproduced the main features seen in the Fig. 4 numerical map, in particular, the cascade to chaos in the larger period-doubling region containing ω = 1. However, there was some variation in the details. For example, the narrower left-side period-doubling cascade in Fig. 8 located near ω = 0.5 shows a small island of the third PD. Some diodes only had the second doubling to period-4 LC before reverse doubling back to period-2 LC. And 3 of the 15 diodes missed this narrow second doubling regime entirely. These 3 diodes all had relatively small reverse-recovery times out of the 15 diodes.

Simulations and circuit measurements were compared for the driven nonlinear RLC-diode circuit where the nonlinearity is due to a rectifier diode in parallel with a fixed capacitor. Three nonlinear characteristics of the diode were included in the simulations: the junction capacitance, the diffusion capacitance, and the DC current/voltage relation. Bifurcation analysis and Lyapunov exponents were used to predict the limit cycles, their stability, and regions of chaos throughout the 2-dim parameter space of the driving source’s frequency and amplitude. Circuit measurements were made in which the function generator’s frequency and amplitude were varied so as to collect data in the same 2-dim parameter space as the simulations. The results are shown in the regime maps, Figs. 4 and 8.

These two regime maps show significant agreement of simulation and circuit measurement over the 2-dim parameter space. The basic structure of the maps is the same except for the sharp dip of the first PD border at twice the resonant frequency. Except for this dip, both maps show that for a low source amplitude, the main LC (frequency the same as the source) covers the entire frequency range, down to the lowest investigated frequency (about 0.2). The bottom panels of Fig. 9 show the time-series for simulation and measurement of this LC for a low source amplitude at ω = 1. At a low source amplitude, the capacitance is essentially constant: the sum of the explicit capacitor, the stray capacitance, and the diode’s junction capacitance. As expected for a constant capacitance with source at resonance, the capacitor–diode voltage amplitude (orange) is much larger than the source (blue) and is nearly sinusoidal. Also, varying the source frequency at a low source amplitude produces the resonance peak expected for an underdamped driven RLC circuit, as seen for the 25 mV source in Fig. 3.

A region of LPs of the main LC is apparent in both maps beginning at a low source amplitude near resonance ( ω = 1) and extending to a lower frequency and a higher source amplitude. These LPs are due to the hysteresis of the ω-continuation seen in Fig. 2. As the source amplitude increases, the resonance peak distorts, creating hysteresis. This distortion is apparent in Fig. 3 where the left side of the peak for the larger source amplitude is getting very steep.

The regime maps in Figs. 4 and 8 show similar period-doubling cascades to chaos. Time series related to the period-doublings are seen in Fig. 9 where the transition from the bottom panel to the middle panel corresponds to crossing the first PD line in the regime maps as the source voltage increases at ω = 1. The resulting period-2 LC for simulations and measurements has the same patterns. The transition from the middle panel to the top panel corresponds to crossing the second PD line, which results in period-4 LC. Again, agreement on the LC shape is good between the left and right panels. Chaos inside the period-doubling cascades of the maps is predicted by the regions of positive maximal LE in Fig. 6. Locations of PD lines in Fig. 4 agree with the locations indicated in Fig. 6. The phase plots for simulations and measurements in Figs. 7 and 10(b) look quite similar and are also consistent with chaos inside the period-doubling cascades.

The agreement discussed above between simulations and circuit measurements of the dynamical regimes suggests that the standard diode model used here is robust. It was necessary to take all three nonlinear characteristics of the diode into account in order to achieve this level of agreement. The one major difference between simulations and circuit measurements in Figs. 4 and 8 is the sharp dip of the first PD bifurcation at ω = 2 seen in simulation but not in measurement. This same discrepancy is apparent in earlier comparisons of simulations and measurements made by Carroll and Pecora.11 They used measurements from different diodes to make regime maps showing the onset of the first PD and then compared the maps to calculated first PD regime maps. Only some of the diode circuit measurements showed the PD dip at ω = 2, whereas nearly all the simulations showed the dip. Their result is consistent with the finding here that the measurements in Fig. 8 seem to show an upside-down shoulder at ω = 2, but not a dip, whereas the simulation in Fig. 4 has the dip. This discrepancy of the existence of the PD dip was not a concern in the earlier work and is an interesting observation here to be investigated in future work.

Other early works focused on identifying the nonlinear property of the diode responsible for the period-doubling cascades to chaos. Accounting for both the nonlinear junction capacitance and the diffusion capacitance predicted observed bifurcations.6,8 Diffusion capacitance can account for the diode’s history-dependent reverse-recovery time, thereby predicting period-doubling cascades to chaos.20 Later, a mathematically simpler model obtained from the junction capacitance alone by using a cubic approximation of Eq. (A8) was used to predict nonlinear phenomena seen in the varactor-loaded SRRs.14,15

Recently, simulations inspired by applications of nonlinear SRRs were done, producing 2-dim regime maps using the maximal LE, which predict period-doubling, chaos, hysteresis, and multistability.21 Those simulations take into account the varactor diode’s junction capacitance by using the cubic approximation mentioned above. There were no experimental results for comparison.

The good agreement of prediction and measurement presented in Figs. 4 and 8 required all three nonlinear characteristics of the diode. Equation (A8) was used instead of its cubic approximation so as to not limit the range of parameter space. This level of accurate prediction could be beneficial for designing varactor-loaded SRR use in metamaterials.

The scope of this study does not include an extensive search for all possible dynamical behaviors. For example, there may be regions of coexisting stable states not presented here. Also, small stable periodic windows within the regions of chaos were not explored.3 These finer detail dynamics may be investigated in future work. Here, the intention is to determine how well the standard diode model works by comparing the dominant dynamical behaviors found in numerical simulations and in circuit measurements over a 2-dim parameter space. The dominant dynamical behaviors have the large basin’s of attraction and are, therefore, not strongly dependent on initial conditions and are the states most easily observed in the circuit.

The diode resonator is a well-known circuit demonstrating nonlinear behavior, including period-doubling cascades to chaos. Earlier works have shown comparisons of simulated and measured bifurcation diagrams; however, comparisons of dynamical regime maps are far less common. This work adds to the long story of analysis of the diode resonator by using bifurcation analysis tools to generate more detailed regime maps for comparison with circuit measurements. The standard diode model used in the SPICE circuit simulator gave good results when used in the bifurcation tools and Lyapunov exponent calculations. The location of period-doubling cascades, hysteresis, and chaos in the 2-dim parameter space of drive frequency and amplitude is nearly the same in the simulated and measured maps. Lyapunov exponent plots are consistent with the regime maps. The bifurcation methods and tools used to generate the regime maps may be useful to the emerging applications of nonlinear SRRs. Future work will address the first PD dip discrepancy in the low-amplitude/high-frequency portion of the regime map. In addition, other values of resistive loss (R in circuit, parameter σ in model) will be investigated.

This work was supported by the University of North Carolina Greensboro Research Assignment.

The author has no conflicts to disclose.

Edward H. Hellen: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are available within the article.

The mathematical model for the semiconductor diode has three parts: the quasi-static (DC) current–voltage relation, a junction capacitance important for reverse bias and small forward bias voltages, and a diffusion capacitance important for accounting for the diode’s reverse-recovery time.20 The model’s parameters are the same as those used in SPICE.1 

The well-known “diode equation” gives the DC current–voltage relation. Its parameters are reverse saturation current I S and emission coefficient N. The mathematical model is
I d ( V d ) = I S ( e 40 V d / N 1 ) ,
(A1)
where 40 comes from e / k T = ( 25 mV ) 1 at room temperature.
The model’s parameters for the junction capacitance are the zero-bias junction capacitance C J 0, the contact potential V J, the junction capacitance grading exponent M, and the forward bias coefficient F C. The junction capacitance’s dependence on voltage is given by
C J ( V d ) = C J 0 ( 1 V d V J ) M .
(A2)
Parameter F C is in the range 0–1 and gives an upper limit of F C × V J for V d in order to avoid the singularity at the forward bias value V d = V J.
The diffusion capacitance needs one additional parameter, the forward transit time T T. The diffusion capacitance’s dependence on voltage is given by
C d f ( V d ) = 40 T T × I S N e 40 V d / N .
(A3)

We use parameter values typical for the 1N4003 rectifier diode: I S = 70 pA, N = 1.4, C J 0 = 33 pf, V J = 0.35 V, M = 0.45, F C = 0.5, and T T = 5 μs. We use the 1N4003 diode for circuit measurements because it has a larger junction capacitance than small signal diodes, thereby reducing the effects of stray capacitance. We measured the diode capacitance as a function of bias voltage for comparison with the model’s predictions. The results in Fig. 11 show our data and the models for the junction capacitance and the diffusion capacitance using parameter values given above. The plot shows that junction capacitance dominates for reverse bias and that diffusion capacitance dominates for forward bias (>0.3 V).

FIG. 11.

Junction capacitance and diffusion capacitance vs bias. Measured data values (X) for 1N4003.

FIG. 11.

Junction capacitance and diffusion capacitance vs bias. Measured data values (X) for 1N4003.

Close modal
The dependence of v d on the variables ( q , i ) is required in Eq. (12). The junction capacitance Eq. (A2) is useful because by definition,
C J = d Q d V d = C J 0 d q d v d .
(A4)
C J is integrated to find q ( v d ), which is then inverted to find v d ( q ).14,15 M = 1 and M 1 are treated separately. The case M = 1 with the condition q ( v d = 0 ) = 0 results in
q ( v d ) = ln ( 1 1 v d ) .
(A5)
Rearranging for voltage and using dimensionless variables finds
v d ( q ) = 1 e q .
(A6)
For M 1, integration with the condition q ( v d = 0 ) = 0 gives
q ( v d ) = 1 1 M [ 1 ( 1 v d ) 1 M ] .
(A7)
Rearranging for voltage finds
v d ( q ) = 1 [ 1 ( 1 M ) q ] 1 1 M .
(A8)
Thus, the dimensionless result for v d ( q ) is
v d ( q ) = { 1 ( 1 ( 1 M ) q ) 1 1 M for  M 1 , 1 e q for  M = 1 ,
(A9)
with the restriction ( 1 M ) q < 1. Figure 12 shows plots of Eq. (A9) for various M-values. It is important to point out that if ( 1 M ) q > 1, then Eq. (A9) gives physically unrealistic complex numbers for V d since the exponent 1 / ( 1 M ) will in general not be an integer. These complex voltages are avoided by incorporating parameter FC. If FC is not included, then complex voltages are avoided only in carefully selected cases, which give integer exponents [like M = 0.8 giving 1 / ( 1 M ) = 5] because a negative number raised to an integral power remains real. The use of these special integral exponent cases requires caution since the imaginary parts of complex number results remain zero for ( 1 M ) q > 1, thereby masking the nonphysical conditions.
FIG. 12.

Dimensionless diode voltage as a function of junction charge for different M-values.

FIG. 12.

Dimensionless diode voltage as a function of junction charge for different M-values.

Close modal
The unphysical singularity of the junction capacitance at V d = V J is avoided by letting C j remain constant for V d > F C × V J, retaining the value C J 0 / ( 1 F C ) M it has at V d = F C × V J. The q-value where this transition occurs is found using v d = F C in Eq. (A7),
q t = 1 ( 1 F C ) ( 1 M ) 1 M .
(A10)
The constant junction capacitance for v d > F C makes the voltage–charge relationship linear. Therefore, Eq. (A8) is modified to maintain a constant slope for q > q t,
v d ( q ) = { 1 ( 1 ( 1 M ) q ) 1 1 M for  q < q t , F C + ( 1 F C ) M ( q q t ) for  q > q t .
(A11)
Equation (3) shows that d v d ( q ) d q is required. Using Eq. (A11),
d v d ( q ) d q = { ( 1 ( 1 M ) q ) M 1 M for  q < q t , ( 1 F C ) M for  q > q t .
(A12)
The dimensionless DC current–voltage relation is
i D ( q ) = i S ( e 40 V j v d ( q ) / N 1 ) ,
(A13)
where Eq. (4) is used to obtain i S, the diode’s dimensionless reverse saturation current.
1.
W. D.
Roehr
,
Rectifier Applications Handbook
(
Motorola, Inc.
,
1993
).
2.
P. S.
Linsay
, “
Period doubling and chaotic behavior in a driven anharmonic oscillator
,”
Phys. Rev. Lett.
47
,
1349
1352
(
1981
).
3.
J.
Testa
,
J.
Pérez
, and
C.
Jeffries
, “
Evidence for universal chaotic behavior of a driven nonlinear oscillator
,”
Phys. Rev. Lett.
48
,
714
717
(
1982
).
4.
M.
Feigenbaum
, “
Quantitative universality for a class of nonlinear transformations
,”
J. Stat. Phys.
19
,
25
52
(
1978
).
5.
M.
Feigenbaum
, “
The universal metric properties of nonlinear transformations
,”
J. Stat. Phys.
21
,
669
706
(
1979
).
6.
A.
Azzouz
,
R.
Duhr
, and
M.
Hasler
, “
Transition to chaos in a simple nonlinear circuit driven by a sinusoidal voltage source
,”
IEEE Trans. Circuits Syst.
30
,
913
914
(
1983
).
7.
C.
Jeffries
and
J.
Perez
, “
Observation of a Pomeau-Manneville intermittent route to chaos in a nonlinear oscillator
,”
Phys. Rev. A
26
,
2117
2122
(
1982
).
8.
S. D.
Brorson
,
D.
Dewey
, and
P. S.
Linsay
, “
Self-replicating attractor of a driven semiconductor oscillator
,”
Phys. Rev. A
28
,
1201
1203
(
1983
).
9.
T.
Klinker
,
W.
Meyer-Ilse
, and
W.
Lauterborn
, “
Period doubling and chaotic behavior in a driven Toda oscillator
,”
Phys. Lett. A
101
,
371
375
(
1984
).
10.
S.
Tanaka
,
T.
Matsumoto
, and
L.
Chua
, “
Bifurcation scenario in a driven R- L-diode circuit
,”
Physica D
28
,
317
344
(
1987
).
11.
T. L.
Carroll
and
L. M.
Pecora
, “
Parameter ranges for the onset of period doubling in the diode resonator
,”
Phys. Rev. E
66
,
046219
(
2002
).
12.
J.
Pendry
,
A.
Holden
,
D.
Robbins
, and
W.
Stewart
, “
Magnetism from conductors and enhanced nonlinear phenomena
,”
IEEE Trans. Microwave Theory Tech.
47
,
2075
2084
(
1999
).
13.
I. V.
Shadrivov
,
S. K.
Morrison
, and
Y. S.
Kivshar
, “
Tunable split-ring resonators for nonlinear negative-index metamaterials
,”
Opt. Express
14
,
9344
9349
(
2006
).
14.
B.
Wang
,
J.
Zhou
,
T.
Koschny
, and
C. M.
Soukoulis
, “
Nonlinear properties of split-ring resonators
,”
Opt. Express
16
,
16058
16063
(
2008
).
15.
N.
Lazarides
,
V.
Paltoglou
, and
G. P.
Tsironis
, “
Nonlinear magnetoinductive transmission lines
,”
Int. J. Bifurcation Chaos
21
,
2147
2159
(
2011
).
16.
A.
Rose
,
D.
Huang
, and
D. R.
Smith
, “
Controlling the second harmonic in a phase-matched negative-index metamaterial
,”
Phys. Rev. Lett.
107
,
063902
(
2011
).
17.
B.
Ermentrout
, Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students, Software, Environments and Tools (Book 14) (SIAM, 2002).
18.
E. J.
Doedel
,
T. F.
Fairgrieve
,
B.
Sandstede
,
A. R.
Champneys
,
Y. A.
Kuznetsov
, and
X.
Wang
,
Auto-07p: Continuation and Bifurcation Software for Ordinary Differential Equations
(
Concordia University
,
2012
).
19.
A.
Wolf
,
J. B.
Swift
,
H. L.
Swinney
, and
J. A.
Vastano
, “
Determining Lyapunov exponents from a time series
,”
Physica D
16
,
285
317
(
1985
).
20.
R. M.
de Moraes
and
S. M.
Anlage
, “
Unified model and reverse recovery nonlinearities of the driven diode resonator
,”
Phys. Rev. E
68
,
026201
(
2003
).
21.
G. D.
Leutcho
,
L.
Woodward
, and
F.
Blanchard
, “
Nonlinear dynamics of a single-gap terahertz split-ring resonator under electromagnetic radiation
,”
Chaos
33
,
103131
(
2023
).