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.

## I. INTRODUCTION

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 $\u2a85$10s of $\mu $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.

## II. RLC-DIODE EQUATIONS

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.

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.

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$, $FC=0.5$, and $TT=5\mu $s.

^{17,18}The autonomous version of system (5) is then

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

## III. ANALYSIS METHOD

Bifurcation analysis of Eq. (12) was done using XPPAUT^{17} 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 $\omega $-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 $\omega $. 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.

The two LPs in Fig. 2 create a region of hysteresis between their locations at $\omega =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 $\omega $ 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 $\omega $– $\mu $ 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 $(\omega ,\mu )$ parameter space.

$\omega $-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 $\omega $– $\mu $ 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).

## IV. CIRCUIT

Circuits were constructed based on Fig. 1. The inductor was $L=10$ mH with $23\Omega $ 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 $\mu $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\Omega $ making the total series resistance $R=123\Omega $. Using $123\Omega $ 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 $\sigma =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\xd725 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\xd7100 mV=0.5$ V.

## V. RESULTS

### A. Numerical bifurcation analysis

Figure 4 is the numerical dynamical regime map showing the limit cycles and chaotic oscillations of the RLC-diode circuit in the 2-dim $\omega $– $\mu $ parameter space. The map consists of LP and PD lines. Traversing the map at $\mu =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 $\omega $-continuation of Fig. 2.

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 $\mu =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.

### B. Lyapunov exponents

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 $\mu =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 $\omega =0.5$, and a wider region around $\omega =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 $\mu =4$ in Fig. 4.

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 $\omega =1.25,\mu =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.

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 $\omega =0.66$ for $\mu =4$. The third PD is very close to the second, suggesting that the cascade to chaos should also occur close to $\omega =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 $\omega =0.61$.

### C. Circuit measurements

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\Omega $ ( $100\Omega $ resistor added to inductor’s $23\Omega $) 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.

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 $\omega =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.

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 $\omega =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.

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 $\omega =1$. However, there was some variation in the details. For example, the narrower left-side period-doubling cascade in Fig. 8 located near $\omega =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.

## VI. DISCUSSION

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 $\omega =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 ( $\omega =1$) and extending to a lower frequency and a higher source amplitude. These LPs are due to the hysteresis of the $\omega $-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 $\omega =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 $\omega =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 $\omega =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 $\omega =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.

## VII. CONCLUSION

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 $\sigma $ in model) will be investigated.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX: DIODE MODEL AND DERIVATION OF $ v d(q)$

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}

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$, $FC=0.5$, and $TT=5\mu $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).

^{14,15}$M=1$ and $M\u22601$ are treated separately. The case $M=1$ with the condition $q( v d=0)=0$ results in

## REFERENCES

*Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students*, Software, Environments and Tools (Book 14) (SIAM, 2002).

*Auto-07p: Continuation and Bifurcation Software for Ordinary Differential Equations*