Symmetry-breaking transitions in networks of nonlinear circuit elements

We investigate a nonlinear circuit consisting of N tunnel diodes in series, which shows close similarities to a semiconductor superlattice or to a neural network. Each tunnel diode is modeled by a three-variable FitzHugh-Nagumo-like system. The tunnel diodes are coupled globally through a load resistor. We find complex bifurcation scenarios with symmetry-breaking transitions that generate multiple fixed points off the synchronization manifold. We show that multiply degenerate zero-eigenvalue bifurcations occur, which lead to multistable current branches, and that these bifurcations are also degenerate with a Hopf bifurcation. These predicted scenarios of multiple branches and degenerate bifurcations are also found experimentally.


Introduction
Nonlinear circuit elements, in particular semiconductor nanostructures, are known to exhibit a wide range of complex spatio-temporal patterns. Among these structures, the superlattice [1, 2, 3, 4, 5, 6, 7, 8] is a prominent example. It consists of alternating layers of two semiconductor materials with different band gaps. As a consequence this leads to an energy band scheme with a periodic sequence of potential barriers and quantum wells. The electrons are localized in the quantum wells if the potential barriers are sufficiently thick. This gives rise to sequential resonant tunneling of electrons between the individual quantum wells. Applying a dc bias, various nonlinear charge transport phenomena may occur, including negative differential conductance (NDC). The nonlinearity may induce spatio-temporal patterns, examples being stationary or traveling high and low field domains. Stationary domains are manifested as multistable branches in the currentvoltage (I-V ) characteristics, where each branch corresponds to the localization of the domain wall in a certain quantum well of the superlattice. In this work we focus on current branches, which form saw-tooth like patterns in the I-V characteristic when the voltage is swept up or down [9,10,11,12,13,14,15,16].
We aim to understand the phenomenon of current branches occurring in superlattices and other semiconductor nanostructures by constructing a breadboard equivalent to the basic electronic transport characteristics in these devices. In fact, we will show that the electronic transport in superlattices, which leads to current branches, can be approximated by an equivalent circuit consisting of tunnel diodes connected in series. Tunnel diodes are nonlinear circuit elements, with similar N -shaped currentvoltage characteristics as superlattices and double barrier resonant tunneling diodes (DBRT ) [17,18,3,19,20,21,22], which can also show a regime, where the current I decreases with increasing voltage V . This is depicted in Fig. 1, where the NDC regime of the schematic I-V curve is shaded for illustration. If the tunnel diode is operated in a circuit with bias voltage V 0 and load resistor R, the load line (dashed black) is given by I = (V 0 − V )/R by Kirchhoff's laws.
In this work we investigate a chain of tunnel diodes connected in series which results in a global (mean field) coupling. The circuit equations, which we derive in Sec. 2, model the NDC elements by a three-variable extension of the simplified FitzHugh-Nagumo (FHN) system. The model reduces to the latter for a certain choice of circuit parameters. Already in 1962, Nagumo et al. pointed out that a model of a single neuron introduced by FitzHugh [23] can be approximated by an electrical circuit containing a tunnel diode [24]. Therefore, our results do not only have important applications to electronic transport in NDC devices such as superlattices, but also to the dynamics of neurosystems [25,26].
The aim of this work is to study the bifurcations in a network of tunnel diodes, both theoretically and experimentally. At the same time this will offer insight into the dynamics of neural networks and electronic transport in superlattices. We find Hopf bifurcations degenerate with zero-eigenvalue bifurcations that are pitchfork bifurcations for two coupled tunnel diodes, and pitchfork or transcritical bifurcations for larger numbers of coupled tunnel diodes. These degenerate bifurcations represent symmetrybreaking nonequilibrium phase transitions [42], where the symmetric, synchronized state becomes unstable and two stable skew-symmetric states arise. These symmetry-breaking bifurcations result in current branches in the serial array of tunnel diodes: by increasing the bias voltage, each tunnel diode subsequently passes through the NDC regime, where oscillations arise. Hopf-pitchfork bifurcations, which we find in two coupled tunnel diode elements, were recently observed in a spatio-temporal FHN-system [43] and in a pacemaker system [44] and were investigated theoretically in [45,46,47,48]. The reason for the simultaneous occurrence of Hopf and pitchfork bifurcations in our system of coupled tunnel diodes is the coexistence of a reflection symmetry and identical Hopf bifurcations in each of the single 3-variable tunnel diode systems.
Pitchfork bifurcations usually arise in systems with Z 2 -symmetry and are generic in these systems. Introducing a heterogeneity breaks the symmetry and leads to an imperfect pitchfork bifurcation, i.e., in systems without a Z 2 symmetry pitchfork bifurcations are of codimension-two. The Hopf-pitchfork bifurcation can be unfolded by separating Hopf and pitchfork bifurcation or by lifting the degeneracy of the pitchfork bifurcation, which we both demonstrate for our system below. In this sense the Hopfpitchfork bifurcation is a codimension-three bifurcation. Our results from analytical considerations and numerical simulations are verified by measurements of tunnel diode circuits.
This paper is organized as follows. In Sec. 2, the circuit and the model will be introduced and compared with electronic transport models for superlattices and with the FHN system. In Sec. 3, a single element will be investigated analytically and numerically. Section 4 deals with a circuit of two tunnel diodes in series, and Sec. 5 treats an arbitrary number of tunnel diodes, where we find a generalization of the Hopfpitchfork bifurcations found in Sec. 4 depending on the number of coupled tunnel diode elements. In Sec. 6, measurements are compared to numerical simulations of the tunnel diode circuits. Finally, conclusions are presented in Sec. 7.

The circuit model
Figure 2(a) shows the circuit of N tunnel diodes connected in series. A single element, which contains the internal device capacitance C d and inductance L d of the tunnel diode [49,3] as well as an external capacitor C, is shown in Fig. 2(b). The superscript n labels the n-th circuit element. The current f (V d ) represents the ideal tunnel diode I-V characteristic, shown as solid black curve in Fig. 1, by neglecting any device capacitance and inductance. The I-V curve f (V d ) of an ideal tunnel diode can be approximated by a cubic polynomial function [24] f where V d and f (V d ) are the voltage across and the current through the tunnel diode, respectively. The parameters V d 0 and I 0 are the voltage and the current at the inflection point of the cubic function, respectively, K is the distance of V d 0 to both extrema and −1/ρ is the slope at V d 0 , as shown in Fig. 1.
Coupling N tunnel diode elements in series as shown in Fig. 2(a) yields the following circuit equations using Kirchhoff's laws: where i = 1 , . . . , N . The variables I and I L dİ

Tunnel diodes as a superlattice model
Let us discuss the similarities of superlattices and tunnel diodes by deriving the regime of validity of Eqs. (3a)-(3c) as a model for the carrier dynamics in a superlattice following Refs. [6,13,15]. A superlattice consists of a sequence of quantum wells i = 1, . . . , N and can be described by a system of differential equations for the electrical field strength F i at the ith potential barrier. Here, is the permittivity and the term dF i /dt describes the displacement current density. The conduction current density J i→i+1 (F i ) between well i and well i + 1 is a tunnel current and can roughly be approximated by a cubic polynomial as in Eq. (1). The sum of the displacement current density and the conduction current density gives the total current density J, which is independent of the well index.
In the tunnel diode model the displacement current corresponds to the current through the capacitor C and the conduction current corresponds to the current through the inductor L d . The equation resulting from Kirchhoff's laws is Eq. (2b). In the approximation of vanishing device inductance (L d = 0) and device capacitance (C d = 0) the current I (i) L is given by the current-voltage characteristic of the tunnel diode f (V (i) ) (see Eqs. (3b) and (3c)) and Eqs. (3a)-(3c) reduce to: The dynamical equation (4) describing a superlattice is formally identical to the dynamical system (5) describing a reduced tunnel diode circuit. A more intuitive approach to explain the similarity between both dynamical systems is to model each potential barrier in a superlattice by a tunnel diode with a parallel capacitor, which describes the charge accumulation in the quantum well, resulting in building up the field in the barrier. The quantum wells are in this case reduced to wires connecting the individual tunnel diode-capacitor circuits to each other. Equation (5) is only an approximation of the carrier dynamics in superlattices since it uses only a very rough approximation of the conduction current density J i→i+1 (F i ). The conduction current density introduces an additional local coupling besides the global coupling through the circuit in Eqs. (4) and (5), since J i→i+1 (F i ) depends also upon the carrier densities n i , n i+1 in the neighboring quantum wells, which are in turn coupled to the fields F i by Gauss' law [3]. This additional nearest-neighbor coupling cannot be observed in the tunnel diode model. Therefore, spatio-temporal patterns like high and low field domains cannot be modeled by the tunnel diode circuit.

Rescaling the model
In the following, a scaling is introduced to simplify the circuit model (3a)-(3c) for all numerical and analytical investigations in this paper. The scaling combines circuit parameters and eliminates the physical dimensions from the dynamical variables so that a set of only four independent dimensionless parameters remains for the analysis. Using this scaling, the analogy of the tunnel diode with the FitzHugh-Nagumo model, which is widely used to model neuronal dynamics, becomes more apparent.
Considering identical tunnel diode elements, we apply the following scaling. The dimensionless variables and parameters are given by With these new variables, the system (3a)-(3c) simplifies to Only the four dimensionless parameters a (N ) , d, , and γ remain in the scaled equations, where and γ are timescale parameters, given by the reactive circuit components. The parameter d is given by the ratio of slopes of the load line and the tunnel characteristics at the inflection point (see Fig. 1) and will be limited to 0 < d < 1, such that there is a single intersection of the load line and the tunnel diode I-V curve, resulting in a single fixed point. The parameter a (N ) depends on the number N of elements and contains the bias voltage V 0 . The calculations in this work are performed with dimensionless variables, but for better comparison with experimental data all results are presented in terms of the circuit variables V L and V (i) , corresponding to x i , y i and z i , respectively.
In Tab. 1 the circuit parameters of the measurements and simulations are listed. The simulation parameters are obtained by using the best possible fit of the measured data, by considering the fabrication variance of the devices. This leads to the small difference in parameters for C and R. In Tab. 2 the physical units of the dimensionless variables derived from the parameters are given. A tunnel diode modeled by Eqs. (7a)-(7c) reduces to the FHN model [23,24,25] for R = 0, which leads to γ = 0 and d = 0. In the case of a single element (N = 1), the = 19.23 τ = tL/ρ = (0.116t) ns γ = 23.18 Table 2. Scaled parameters and variables derived from Tab. 1 using the scaling (6).
parameter a (1) is then equal to z. The resulting system reads These equations are of the same form as the simplified FHN system, widely used in neurodynamics, e.g., [50,51,52,53,54,55]. Note that the values of the parameters and γ we have used to describe the experiment differ from the ones commonly used in the FHN system. We stress, however, that the location of the fixed points is independent of these parameters. Furthermore, for other circuit parameters small and γ can be realized, resembling the FHN model.

A single tunnel diode
First, we perform a bifurcation analysis of a single tunnel diode. The dynamics of a single diode resembles the synchronized dynamics of N tunnel diodes in series, however with a different parameter a (N ) , see Eqs. (6). The fixed point (x * , y * , z * ) of Eqs. (7a)-(7c) for N = 1 is given by a third order polynomial for the z component: and x * = −z * , y * = (a − z * )/d. The fixed point depends on the parameter a (1) , which includes the applied voltage V 0 , and on the ratio d = R/ρ (see Eqs. (6)). To calculate the stability of the fixed point, we consider the eigenvalues of the Jacobian   For each curve, i.e., each value of d, the fixed point is unstable in a specific interval of V 0 . Using the eigenvalues of the Jacobian matrix and numerical bifurcation analysis we will show that the fixed point changes its stability at the boundaries of this interval in subcritical Hopf bifurcations.
Figures 3(c) and (d) show a bifurcation diagram at the first subcritical Hopf bifurcation (low V 0 ) for d = 0.57. The plots were obtained using the continuation tool Auto [56]. The light blue (stable) and dashed dark blue (unstable) lines show the maxima and minima of V (1) (c) and I (1) L (d) of the periodic orbit. For better visualization an ellipse is plotted for each branch of stable or unstable periodic orbits. The position where the subcritical Hopf bifurcation (HB) occurs is denoted by the blue dot at V 0 = 0.2698 V. At V 0 = 0.2640 V the unstable and stable limit cycles collide in a saddle-node bifurcation of limit cycles. An illustration of the subcritical Hopf bifurcation in three dimensional phase space can be found in movie 1 in the supplementary material.
For very small values of the resistance R, i.e., in the limit of the simplified FHN system, the Hopf bifurcations occur at the extrema of the current-voltage curve (blue curve in Fig. 3(b)). This shows that the tunnel diode system gives a smooth approach to the simplified FHN system.
In the following sections we will show that the synchronized state of N tunnel diodes shows similarities to a single tunnel diode but with a stretching factor N as introduced in parameter a (N ) in Eqs. 6.

Two tunnel diodes in series
In this section we perform a bifurcation analysis of two tunnel diodes connected in series. We see that a codimension-three Hopf-pitchfork bifurcation leads to a symmetrybreaking transition. Introducing heterogeneity to the system of tunnel diodes lifts the degeneracy and unfolds the Hopf-pitchfork bifurcation.

Identical tunnel diodes
Similar to the analysis of a single element in Sec. 3, we can calculate the fixed points in the case of N = 2 in Eqs. (7a)-(7c). We find two coupled third order polynomial equations for z * 1 and z * 2 of the fixed point ( Again both equations depend on a (2) (V 0 ) and on the ratio d = R/ρ. The Jacobian at the fixed points is given by: The stability of the fixed points depends on the x component of both fixed points, as well as the timescale parameters and γ and parameter d.  The point symmetry discussed above for a single element (refer to Sec. 3) remains valid for the system of two tunnel diodes. But the symmetry point is shifted to (V 0 , V ) = (0.594V, 0.267V) due to the rescaled value of a (N ) .
Another symmetry in the two-diode system is a Z 2 symmetry corresponding to the interchange of the dynamical variables of tunnel diode one and two, i.e., the Eqs. (11a) and (11b) and (7a)-(7c) are invariant under the transformation (x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ) → (x 2 , y 2 , z 2 , x 1 , y 1 , z 1 ) (13) and the phase space is symmetric with respect to reflections at the synchronization manifold. § The synchronization manifold (x 1 , y 1 , z 1 ) = (x 2 , y 2 , z 2 ) shows identical behavior to the dynamics of a single tunnel diode only with a rescaled parameter a (N ) as introduced in Eqs. (6). All fixed points represented by the middle branch in Fig. 4 are located within the synchronization manifold, since both tunnel diodes voltages V (1) and V (2) are the same.
The pitchfork bifurcations at V 0 = 0.334V and V 0 = 0.855V for d = 0.01 and at V 0 = 0.419V and V 0 = 0.619V for d = 0.99 mark a symmetry-breaking transition: the synchronized state becomes unstable. Note that a Hopf bifurcation coincides with the pitchfork bifurcation, thus the two emerging fixed point branches are unstable and the system will actually follow one of two stable periodic orbits, as we will show below.
The pitchfork bifurcation can be super-or subcritical depending on the parameter d as shown in Fig. 4. The change from a subcritical pitchfork bifurcation to a supercritical pitchfork bifurcation occurs at d = 0.34 for the parameter set in Fig. 4.
For all values of d between zero and one, a Hopf bifurcation occurs simultaneously with the pitchfork bifurcations, i.e., a complex conjugate pair of eigenvalues as well as a real eigenvalue cross the imaginary axis in the complex plane simultaneously. To find the location of this Hopf-pitchfork bifurcation further, we proceed as follows: The characteristic polynomial χ(λ) is of 6th order and thus difficult to solve directly. However, below the Hopf-pitchfork bifurcation there is only one symmetric fixed point, i.e., x ≡ x 1 = x 2 . Furthermore, at the pitchfork bifurcation we have χ(0) ! = 0 which gives (after calculating the characteristic polynomial)

This equation has four solutions for
Since we consider 0 < d < 1, the latter two solutions are imaginary and thus spurious and we have found the coordinate values x = ±1 at the Hopf-pitchfork bifurcations. The fixed point equations (11a) and (11b) reduce for Inserting now x = ±1 we find the location of the Hopf-pitchfork bifurcation In Appendix A.1, we analyze the characteristic polynomial in the general case of N tunnel diodes. § As introduced in [57] the synchronization manifold is a hyperplane in phase space of a system, where the synchronized dynamics take place. Figure 5 shows a bifurcation diagram at the first subcritical pitchfork bifurcation (low V 0 ) for d = 0.57 obtained with the continuation tool Auto. Again, solid black lines represent stable fixed points, while dashed black lines represent unstable fixed points. The bifurcation diagram shows possible solutions for one of the tunnel diodes only. If one of the tunnel diodes is located on the middle (symmetric) fixed point, the other tunnel diode is also in this fixed point. If one of the tunnel diodes is located on the upper branch of the fixed points, the other tunnel diode is on the lower branch and vice versa.
The green, orange, and blue lines show the maxima and minima of periodic orbits, where solid lines and light color shading denote stable periodic orbits, and dashed lines and dark color shading denote unstable periodic orbits. For better visualization, an ellipse is plotted for each branch of stable or unstable periodic orbits. We will discuss the differences between the green, orange, and blue orbits below. In Fig. 5 the Hopf-pitchfork bifurcation occurs at V 0 = 0.381V. Additionally, at each fixed point branch of the pitchfork bifurcation another Hopf bifurcation occurs for larger values of V 0 : at the two outer branches at V 0 = 0.394V and at the inner branch at V 0 = 0.492V.
The phase space symmetry Eq. (13) implies that there can be three different types of periodic orbits, which we call synchronized orbits, anti-phase orbits, and skew orbits, respectively. The synchronized orbits (shown in blue) lie in the synchronization manifold and thus obey (x 1 (t), y 1 (t), z 1 (t)) = (x 2 (t), y 2 (t), z 2 (t)).
The anti-phase orbits (orange) wind around the 3-dimensional synchronization manifold in the 6-dimensional phase space and obey the strict anti-phase relation where T is the period of the orbit. Finally, the skew orbits (shown in green) are neither anti-phase nor synchronized and come in pairs, due to the phase-space symmetry Eq. (13). They are similar to the outer pitchfork branches for fixed points, i.e., one of the diode voltages oscillates around a high voltage and the other around a low voltage state.
At V 0 = 0.434V the anti-phase orbit surrounding the inner unstable fixed point becomes stable in a pitchfork bifurcation of periodic orbits (PPO) and generates a pair of unstable skew orbits. As is apparent from this bifurcation the skew orbits (at least close to the pitchfork bifurcation) show a generalized anti-phase dynamics, i.e., they do not have a strict anti-phase relation such as the anti-phase orbits, but when there is a maximum in the time-series of the high voltage diode, there is a minimum in the time-series of the low voltage diode and vice versa. In our numerical simulations we see this behavior for all skew orbits and all parameters. We illustrate this behavior in movie 3 which can be found in the supplementary material.
The periodic orbits surrounding the middle fixed point branch in Fig. 5 are strict anti-phase orbits and both tunnel diode voltages oscillate in complete anti-phase. The supplementary movie 4 illustrates this behavior. Figure 6 is a schematic three-dimensional projection of the six-dimensional phase space showing the codimension-three Hopf-pitchfork bifurcation, with values of V 0 as indicated in Fig. 5 as blue lines. The arrows illustrate the direction of the stable and unstable manifolds and thus the stability of the limit cycles and fixed points. An illustration of the complete bifurcation scenario shown in Fig. 5 can be found in the supplementary material in movie 2. The investigation of the Hopf-pitchfork bifurcation by a schematic projection follows the same spirit as the figures in [58], where the stabilization of complex spatio-temporal dynamics near a subcritical Hopf bifurcation is investigated by time-delayed feedback.
In Fig. 6(a) a voltage of V 0 = 0.29V is applied. Only one stable fixed point exists. At V 0 = 0.34V, as in Fig. 6(b), a stable and an unstable skew orbit are born in a saddle-node (fold) bifurcation on each side of the middle stable fixed point. The stable skew orbits move away from the stable fixed point while the unstable skew orbits move towards the fixed point with increasing V 0 .
In Fig. 6(c) at V 0 = 0.38V, which is just below the value of the Hopf-pitchfork bifurcation, two fixed points are born in a fold bifurcation on each side of the middle stable fixed point: one unstable fixed point and a saddle point, which is stable in the direction of the pitchfork bifurcation but unstable in the transversal direction. Also, a fold bifurcation of anti-phase orbits takes place around the stable middle fixed point. Both limit cycles are unstable in the direction of the pitchfork bifurcation.
At the Hopf-pitchfork bifurcation two unstable skew orbits and one of the unstable anti-phase orbits as well as both unstable fixed points close to the center fixed point coalesce with the latter. The stable fixed point in the center becomes unstable. Figure 6(d) illustrates this schematically. The outer two stable skew orbits and two unstable fixed points as well as the unstable anti-phase orbit remain unchanged.
In Fig. 6(e), at V 0 = 0.42V, the two outer, unstable fixed points have become stable in a subcritical Hopf bifurcation by giving birth to two unstable skew orbits. These orbits are moving towards the stable skew orbits which themselves are moving away from the unstable fixed point in the center.
Finally, in Fig. 6(f) at V 0 = 0.45V the two outer unstable skew orbits collide with the two stable skew orbits in a saddle-node bifurcation and the limit cycle around the centered fixed point becomes stable in a PPO bifurcation and gives rise to two unstable skew orbits.

Introducing heterogeneity
So far both tunnel diodes were assumed to be identical as assumed in Eqs. (7a)-(7c). This assumption is not justified in experimental situations (see Sec. 6), since it is impossible to find two perfectly identical electronic devices. In what follows we investigate nonidentical devices by introducing heterogeneity to the parameters of the subsystem. The heterogeneity may be introduced as discrepancy in the parameter a (2) or as discrepancy in one of the time scale parameters or γ in Eqs. (7a)-(7c): In Fig. 7(a) the time scale parameter is chosen as 1 = 19.23 for the first subsystem and 2 = 21.15 for the second subsystem.
The time-scale parameters do not influence the position of the fixed points (see Eqs. (7a)-(7c)) and thus the curves in Fig. 7(a) and Fig. 4 are identical. Only the stability of the fixed points changes and the Hopf bifurcation point shifts. This separates the Hopf and the pitchfork bifurcation from each other. A similar result can be obtained by introducing heterogeneity in the parameter γ. Fig. 7(b) shows this unfolding of the Hopf-pitchfork bifurcation due to a difference ∆ = 2 − 1 in the parameter . The dashed line corresponds to a pitchfork bifurcation while the red, green and blue lines mark the Hopf bifurcation curves in the (V 0 , ∆ )plane. Each color represents a different value of the difference in the timescale parameter ∆γ = γ 2 − γ 1 . In each case a Hopf-pitchfork bifurcation can be obtained by using the correct set of parameter values.
A mismatch in a (2) , on the other hand, does influence the location of the fixed points and in fact breaks up the pitchfork bifurcation and leads to an imperfect pitchfork bifurcation as discussed in [60]. The pitchfork bifurcation separates into a stable branch and a saddle-node bifurcation. Figure 7(c) shows this imperfect bifurcation. The parameter a (2) was previously set to a (2) = 7.6570V 0 − 4.5468 for both elements. This The value is calculated by introducing a difference of 10% in the tunnel diode parameter C d , which is the typical variance for general purpose tunnel diodes as denoted in the data sheet in [59].  In the next sections we give an outlook into the dynamics of N coupled tunnel diodes by performing a bifurcation analysis for a system of four tunnel diodes in detail, and comparing the results to electron transport models in superlattices. ¶ The value is calculated by introducing a difference of 1% in the tunnel diode parameter v d0 , which is the typical variance for high quality tunnel diodes.

Arbitrary number of tunnel diodes in series
In this section we perform a bifurcation analysis of four identical tunnel diodes connected in series and aim to understand a general serial array of N tunnel diodes. Due to the high symmetry the previously described Hopf-pitchfork bifurcation becomes a more degenerate [61] symmetry-breaking transition and we have to distinguish between N being even or odd. These highly degenerate bifurcations lead to multistability between different fixed point branches and show close similarities to current branches in superlattices. This similarity will be addressed below.
The linear stability analysis can be carried out in the same fashion as for one and two tunnel diode elements. The general system for the z * n (n = 1, . . . , N ) values of the fixed points (x * 1 , y * 1 , z * 1 ,. . .,x * N , y * N ,z * N ) is given by by using Eqs. (7a)-(7c). System (14) contains N equations with n = 1, . . . , N . The Jacobian at the fixed point is a block matrix with where the blocks C account for coupling contributions, and x * n = −z * n (n = 1, . . . , N ). The N -dimensional system of equations for the fixed points Eq. (14) and the 3N -dimensional Jacobian (15) can be solved numerically for a given value of N . In the following, we present the results for four tunnel diodes in series with the aim to understand the experimental results in the next section. The bifurcation pattern for an arbitrary number of tunnel diodes is then explained based on these results.  Apart from the symmetric branch, where all z-values are equal, there are two other types of branches: Branches which have only two species, i.e., each tunnel diode can be in either of two different states z, are called primary branches. Solutions with three species, on the other hand, are called secondary branches. Primary branches emanate from the symmetry-breaking point and are unstable for N ≥ 3 close to the bifurcation point (even without the additional Hopf bifurcations in our particular case). Primary branches can become stable in two ways: They can bend back at a fold point (green squares in Fig. 8(a)) or they can be stabilized through bifurcations with secondary branches (yellow diamonds in Fig. 8(a)). In our case due to the additional Hopf bifurcations at the symmetry-breaking point, the primary branches have to undergo another Hopf bifurcation in order to become stable again. In the case of stabilization through bifurcations with secondary branches, the additional Hopf bifurcations coincide with the secondary bifurcations (yellow diamonds in Fig. 8). Note that secondary bifurcation points can be found analytically (see Appendix A.2).
Let us now discuss in more detail the branch structure close to the symmetrybreaking transition. In Fig. 8(a) the possible tunnel diode voltages are plotted. In this projection there are two branches of transcritical type visible, i.e., branches which exist below and above the bifurcation point. Additionally there is one branch of pitchfork type with a vertical tangent visible.
On the transcritical branches one tunnel diode is in one state (on the leftmost branches for instance in the high voltage state) and three are in the other state. How many of such solutions exist? Since there can be one of the four tunnel diodes in the single voltage state, there are four such branches. Note that in the projection of Fig. 8 these branches coincide, since they have the same tunnel diode voltages, only in a different order. On the other hand, on the pitchfork branch, two tunnel diodes are in the high voltage and two tunnel diodes are in the low voltage state, respectively. Since there are 4 2 = 6 ways to choose two from four tunnel diodes, there are 6 solutions of this pitchfork type, which again coincide in the projection of Fig. 8.
All of the above discussed solutions meet at the symmetry-breaking bifurcation point. Since a complete bifurcation analysis of all these primary and secondary branches is beyond the scope of this work, we focus on the experimentally observable features as shown in Fig. 8.
The stable branches have the shape of current branches in the current component I (i) L as shown in Fig. 8(b). This means that by increasing the applied voltage V 0 the current I L suddenly drops to a lower current. This drop is related to a jump of the voltage V (i) in one tunnel diode, i.e., from a low voltage branch to a high voltage branch. After the jump the current increases to the same threshold again. This behavior, which forms a sawtooth each time, repeats N times (red arrows in Fig. 8(b)). We note that similar current branches have been reported for a very simple model consisting of a series array of N DC elements without any inductance [63]. The separate voltage jump of each tunnel diode corresponds to the separate passing through the NDC regime and as a result we observe current branches.
The symmetry-breaking transition is highly degenerate (see Appendix A.1) and gives rise to a plethora of branches. Due to the coexistence of Hopf bifurcations each branch in Fig. 8 is surrounded by limit cycles similar to the system of two tunnel diode elements, as shown in Sec. 4. When a current branch becomes unstable with increasing V 0 the system may, instead of switching to the next branch, also switch to a periodic orbit, due to the multistability between periodic orbits and fixed points.
As discussed above for four tunnel diodes, we observe branches of transcritical as well as of pitchfork type. From the investigation of N = 5, 6, 7, we found the following generic feature. For odd N all branches are of transcritical type, while for even N we also observe pitchfork type branches on which there are an equal number of diodes in the upper and in the lower state. It is clear that for odd N the tunnel diodes cannot be divided into two equal numbers and thus no pitchfork type branches are present.
With increasing N the degeneracy of the bifurcation at the emanating point increases due to the higher symmetry of the system and more branches appear. For large N the branch structure becomes very similar to the current branches of superlattices. However, there is one main difference between the two systems: In the tunnel diode network there exists no determined order of switching of tunnel diodes because of the previously discussed transformation invariance. This is in contrast to the electronic transport models of superlattices, where high-and low-field domains cause the current branches. High-and low-field domains are groups of consecutive quantum wells, which are either in a high-field state or in a low-field state. With increasing voltage the highfield domain expands by one superlattice period, i.e., the domain wall moves to the neighboring quantum well and forms the next current branch, and so forth. Thus the quantum wells switch consecutively.
In the next section we compare our analytical and numerical results for circuits containing a single, two, or four tunnel diodes to measurements.

Comparison of experimental and numerical results
In this section we compare the results from Secs. 3, 4 and 5 to experiments. As a setup for the measurements, we use the circuit in Fig. 2 with general purpose devices including tunnel diodes of the type 1N3714, ceramic capacitors with capacitance C = 56 pF±2 pF and ceramic resistors with the resistance R = 47 Ω ± 3 Ω. The parameters of the measurement circuit and the parameters used in the simulations are listed in Tab. 1.
For the measurements with a ramped voltage source (increasing or decreasing input voltage V 0 ), we use the arbitrary-function generator Tektronix model AFG3252. We apply a slow frequency sawtooth-like oscillation of 10 Hz such that we are able to observe stable oscillating or fixed point states at each voltage. To decouple the function generator from the circuit, we use a buffer, circuit based on a general purpose operational amplifier LM741.
We use the digital phosphor oscilloscope Tektronix model DPO7104 in combination with the differential probe TD1000 (bandwidth 1 GHz) for measuring the tunnel diode voltage or the active probe TAP1500 (bandwidth 1 GHz) for measuring the overall current.
In the simulations, a noise source is introduced to approximate noise effects in the measurements. The noise is implemented into the dynamical system by addition of an independent Gaussian white noise term Dξ (i) (t) for each system with D being the dimensionless noise intensity set to D = 0.01. The noise is added to the dynamical equation of the z n variable of each subsystem, which is reasonable as z n is the rescaled voltage V (n) .

A single tunnel diode
The measured and simulated I-V curves and voltages V (1) of a single tunnel diode vs. input voltage V 0 are shown in Fig. 9. The upper panels show the I-V curve obtained from measurements (panel (a)) and simulations (panel (b)). The lower panels show the voltage V (1) from measurements (panel (c)) and simulations (panel (d)). While the I-V curve is plotted for upwards and downwards ramping, the tunnel diode voltage is displayed for upwards ramping only. The shaded area in all four panels shows the range of time-dependent oscillations (range and amplitude of voltage V 0 ) with the lines denoting the time average of the oscillations in this regime.
Movie 5 shows a plot similar to Fig. 9(c), but ramping the bias voltage both upand downwards. This illustrates both the oscillating regime (shaded in Fig. 9) and the differences between up-and downwards ramping caused by hysteresis. This hysteresis is caused by the subcritical Hopf bifurcation as discussed in Sec. 3.
Overall, the simulated results are in good qualitative agreement with the results from the experiments. The general characteristic, which exhibits an increasing current I or voltage V (1) followed by an oscillatory regime and another increasing regime, is identical for measurements and simulations. On the other hand, the oscillatory regime and its amplitude differ. For example, during upwards ramping, the measured oscillations start at V 0 = 0.17 V while in the simulations oscillations occur after V 0 = 0.27 V, or the measured voltage oscillations V (1) reach the negative voltage regime, while the voltage oscillations stay positive in the simulations as indicated in Fig. 9. The different amplitudes of voltage oscillations is due to the non-ideal nature of the experimental voltage source. The differences in the bias voltage V 0 at which the oscillating regime sets on can be explained by the restriction of the third-order polynomial approximation of the I-V curve in the simulations. The parameters, as we fitted the curve, match the I-V curve in the low-and high-V 0 regime, but as a trade-off the location of the extrema of the curve -which correspond to the onset of oscillations -is not perfectly reproduced by the fit as shown in Fig. 1.
The typical I-V characteristic of a single tunnel diode, which is shown in Fig. 9, starts with a single current branch and is followed by a sudden collapse of the tunneling current because the energy levels in a tunnel diode become misaligned. In the case of a load resistor R < ρ such that 0 < d < 1, the collapse of the tunneling current (d) corresponds to oscillations. This means that when the tunnel diode passes through the NDC regime it exhibits oscillations. Simulations and measurements both show hysteresis while passing through the NDC regime. The width of the hysteresis regime is ∆V 0 = 0.02 V ± 0.04 V in experiments and in the simulations ∆V 0 = 0.006 V. These results agree with the previous results gained from the bifurcation analysis, where the overlap between stable periodic orbits and stable fixed points has a width of ∆V 0 = 0.006 V. We attribute the difference in hysteresis ranges in measurement versus simulation to the use of a third order polynomial approximation in the simulations.

Two tunnel diodes
The simulated and measured results of two tunnel diodes in series are shown in Fig. 10. Again the upper panels (a) and (b) show the circuit I-V curve during upwards (red) and downwards (blue) ramping for measurements and simulations, respectively. The lower panels are the tunnel diode voltages V (1) (black) and V (2) (green) during upwards ramping of the input voltage V 0 in measurements (panel (c)) and simulations (panel Movie 6 shows the effect of hysteresis in measurements on two coupled tunnel diodes while ramping the bias voltage V 0 both up-and downwards. The movie is otherwise similar to Fig. 10(c). The drift observable in this movie is not caused by differences between the individual tunnel diodes (cf. Fig. 7(c)), but influenced by the non-ideal (i.e, not quasi-static) ramping in the measurements.
Again, the simulated and the measured results are qualitatively similar. Simulation and measurement show hysteresis in the I-V curve of both tunnel diodes, caused by the subcritical bifurcations, which imply coexisting stable states. For example, during upwards ramping of the applied voltage, the oscillations of the second oscillation regime start in the measurement at V 0 = 0.68 V and ends at V 0 = 0.85 V while during downwards ramping the oscillations occur between V 0 = 0.66 V and V 0 = 0.73 V. In the simulations, the oscillations of the second oscillation regime occur during upwards ramping between V 0 = 0.79 V and V 0 = 0.89 V and during downwards ramping between V 0 = 0.75 V and V 0 = 0.81 V as indicated in Fig. 10.
Two current branches occur in the current-voltage characteristic of two tunnel diodes in series (in panel (a) and (b) of Fig. 10). The tunneling current breaks down after reaching a current threshold (I measure = 2.17 mA ± 0.02 mA, I sim = 2.28 mA) in both tunnel diodes separately at different input voltages. This breakdown corresponds to a pitchfork bifurcation in the bifurcation diagram. Prior to the pitchfork bifurcation, both tunnel diodes are in the identical tunneling state, while after the pitchfork bifurcation one tunnel diode has passed through the NDC regime, while the other tunnel diode is still in the tunneling state. These are the two asymmetric branches of the pitchfork bifurcation. They can be seen in panel (c) and (d) of Fig. 10. The current branches occur simultaneously with oscillations since the NDC regime is unstable, which gives rise to oscillations. These oscillations are related to stable skew orbits as discussed previously.

Four tunnel diodes
The simulated and measured results of four tunnel diodes in series, which are shown in Fig. 11, are similar to the results of two tunnel diodes. Again the upper panels (a) and (b) show the circuit I-V curve during upwards (red) and downwards (blue) voltage ramping for measurements and simulations, respectively. The lower panels are the tunnel diode voltages V (1) (black), V (2) (green), V (3) (yellow) and V (4) (brown) during upwards ramping of the input voltage V 0 in measurements (panel (c)) and simulations (panel (d)).
In the case of four tunnel diodes in series, four current branches occur as shown in Fig. 11. Each tunnel diode switches from a low-voltage state to a high-voltage state, separately and via oscillations. Again, each tunnel diode passes through the NDC regime at different levels of the bias voltage V 0 . Both theoretical and experimental investigations reproduce the described behavior as shown in Fig. 11 also shown in movie 7, which resembles Fig. 11(c) but with up-and downwards ramping. Like in movie 6, a drift is observable, which is caused by the non-ideal (i.e., not quasi-static) ramping in the measurements.
The theoretical model describes the current branches as branches with multiple degeneracy generated at symmetry-breaking transitions. Since a tunnel diode can be at the upper or lower branch in Fig. 8(a), there exists no well-defined order of switching of the tunnel diodes. The switching shown in Fig. 11 is just one possible way determined by slight differences between each tunnel diode. This is different from electronic transport in superlattices, where the current branches are caused by high-and low-field domains and consecutive shifting of the domain wall.

Conclusions
In this paper we have investigated a nonlinear circuit, which shows dynamical scenarios similar to a semiconductor superlattice. As a model for a sequence of tunneling barriers we use a series connection of tunnel diodes, described by a three-variable extended FitzHugh-Nagumo system. We have shown that a single tunnel diode system reduces to the simplified FitzHugh-Nagumo system in the limit case of small load resistance.  Our numerical model predicts two subcritical Hopf bifurcations in good agreement with measurements.
In the case of N tunnel diodes in series, we have shown that the model approximately describes electronic transport in superlattices for certain parameter ranges and exhibits dynamics similar to that of a superlattice. This includes multistable current branches, which result in saw-tooth current-voltage characteristics and hysteresis upon up-and down-sweep of applied voltage. As the cause for the current branches we have found degenerate zero-eigenvalue bifurcations. Each of the branches corresponds to a different diode passing through the NDC regime. The order in which the tunnel diodes pass through the NDC regime is determined by small differences between each of the devices. This is in contrast to semiconductor superlattices, where the domain wall is shifted consecutively through the superlattice.
Additionally, we have observed that each zero-eigenvalue bifurcation occurs simultaneously with a Hopf bifurcation such that limit cycles surround each fixed point branch. We have shown that this degenerate bifurcation is followed by a pitchfork bifurcation of periodic orbits in the case of two tunnel diode elements. We have been able to observe this complex bifurcation scenario in measurements. The switching between different branches of the pitchfork bifurcation is accompanied by oscillations.
The comparison with simulated results is in good qualitative agreement.
The results of our work aim to understand and illustrate the complex bifurcation scenario in a network with high symmetry, i.e., invariance under any permutation of the tunnel diodes. We have shown that two tunnel diodes which individually exhibit Hopf bifurcations can form a Hopf-pitchfork bifurcation when coupled and that this degenerate bifurcation is caused by the reflection symmetry of the system. We have also demonstrated that current branches can even occur in networks without local coupling. A global coupling can already induce zero-eigenvalue bifurcations, which are multiply degenerate and results in multiple current branches. This is in contrast to superlattices where local coupling is always present due to sequential tunneling between neighboring quantum wells.
Our results may also help to understand the complex dynamics in cases where two or more tunnel diodes are used in series, for example in multi-junction solar cells with more than two active solar subcells. and it is immediately clear that such a transition is (N − 1) fold degenerate. Using the explicit form of J 1 and C (see Eqs. (16) and (17)) and solving Eq. (A.5) for x, we obtain i.e., the symmetry-breaking transitions always occur at these values of x (see Fig. 5 and 8). Inserting these x-values into Eq. (14) (z n = −x = ∓1) we can also find the value of the bifurcation parameter a (N ) at the symmetry-breaking bifurcation: To find the overall eigenvalue structure at the symmetry breaking we substitute x = ±1 into Eq. (A.4). This gives eigenvalues Thus there is a Hopf bifurcation coinciding with the zero-eigenvalue bifurcation and both of these are (N − 1) fold degenerate. Note that it is possible to calculate the crossing direction of the eigenvalues in the complex plane by implicit differentiation of χ ⊥ (λ) with respect to a. This shows that the fixed point is transversely stable below the first transition (a < −(N − 2d/3)) and above the second transition (a > + (N − 2d/3)). Of course this shows the stability only close to the transitions.
Let us now discuss the longitudinal stability of the fixed point. The characteristic Eq. (A.3) for the parallel direction evaluated at x = ±1, is given by Since all coefficients are positive, the Routh-Hurwitz criterion, which ensures that all eigenvalues are stable is given by and is always fulfilled (d, N > 0). The synchronized fixed point is thus stable in the parallel direction at the symmetry-breaking transition and the bifurcation is observable for any number N of tunnel diodes. Let us now specifically investigate the case of four tunnel diodes (N = 4). To find the bifurcation parameter at the secondary bifurcations, we need to determine the value of σ at the bifurcation. In order for a secondary bifurcation to occur in the first place, at least two tunnel diodes need to have z = 1 or z = −1. This gives the following possible values for σ σ = 2 · (−1) + 2 · (+2) = 2, (A.12) For other values of N the secondary bifurcation points (a-values) can be found similarly.