Complex Dynamics in a Memcapacitor-Based Circuit

In this paper, a new memcapacitor model and its corresponding circuit emulator are proposed, based on which, a chaotic oscillator is designed and the system dynamic characteristics are investigated, both analytically and experimentally. Extreme multistability and coexisting attractors are observed in this complex system. The basins of attraction, multistability, bifurcations, Lyapunov exponents, and initial-condition-triggered similar bifurcation are analyzed. Finally, the memcapacitor-based chaotic oscillator is realized via circuit implementation with experimental results presented.


Introduction
Memcapacitor is a type of memory device composed of a memristor and meminductor [1], which emerged after the realization of a real memristor prototype [2]. A memcapacitor is actually a nonlinear capacitor with instantaneous responses depending on the internal states and the input signals. Although several possible realizations of the memcapacitor were attempted, e.g., using a micro-electro-mechanical system [2], ionic transport [3], or electronic effect [4], a memcapacitor is not yet available commercially in any form. Therefore, building functional analog memcapacitor models and emulators for computer simulations and laboratory experiments has become urgent, attracting immense interest from both academia and industry.
There have been many papers about the memristor [5][6][7][8][9]. In [6], a novel complex Lorenz system with a flux-controlled memristor is introduced and investigated. An active controller is designed to achieve modified projective synchronization (MPS) based on Lyapunov stability theory. In [7], a new memristor-based hyperchaotic complex Lu system is investigated, where an adaptive controller and a parameter estimator are proposed to realize complex generalized synchronization. Furthermore, the complex dynamics of fractional-order and diode bridge-based memristive circuits are studied in [8] and [9]" respectively. Compared with the reports concerning memristors, references of memcapacitors are relatively fewer. Existing research involves designing a memcapacitor SPICE (Simulation program with integrated circuit emphasis) simulator [10][11][12] and mutator that can transform a memristor to a memcapacitor [13,14]. In [15], the boundary dynamics of a charge-controlled memcapacitor is investigated, where Joglekar's window function is used to describe the nonlinearities of memcapacitor's boundaries. In [16], a mathematical memcapacitor model is introduced and a memcapacitor oscillator is designed, with theoretical and experimental analyses on their basic dynamic characteristics given. In [17], a floating emulator circuit is built, using common off-the-shelf active devices, to mimic the dynamic behaviors of flux-coupled memcapacitors. σ M (t) = q M (t) (1) where v c (t) is the voltage across the memcapacitor and C −1 is the inverse memcapacitance. The symbol q M (t) is the charge going through the memcapacitor at time t, and σ M (t) is the integral of q M (t).
To describe the proposed charge-controlled memcapacitor model, Equation (2) is used, where the memcapacitance depends on the device charge and is changing nonlinearly. Since the meminductor has not been fabricated, in this paper we assume the inverse memcapacitance to be precisely defined as C −1 = a + b|σ M (t)|, so as to obtain: v c (t) = (a + b|σ M (t)|)q M (t) . σ M (t) = q M (t) (2) An emulating circuit is designed, as shown in Figure 1, to realize the above charge-controlled memcapacitor. In this circuit, it is assumed that all components are ideal without losses and the output limitations are set as ±15 V.
The circuit resistors are set as R 1 = R 2 , R 3 = R 4 , and R 5 = R 6 . The operational amplifier U 1 is used to reverse the sign of i, which is the current going into the floating terminal of the memcapacitor. Then, the charge q going through capacitor C 1 can be calculated by integrating i in the time domain. The operational amplifier U 2 constitutes a subtraction circuit, used to extract the voltage across the capacitor C 1 , with output voltage: The operational amplifier U 3 is an integral circuit and its output is: Operational amplifiers U 4 and U 5 construct an absolute-valued circuit with output voltage: Now, based on Equations (3) and (5), the output voltage of the multiplier M 1 can be calculated using: Then, the voltage across this grounded memcapacitor can be written as: Entropy 2019, 21, 188 3 of 14 The operational amplifier U3 is an integral circuit and its output is: where σM(t) is the integral of qM(t).
Operational amplifiers U4 and U5 construct an absolute-valued circuit with output voltage: If the memcapacitor-based emulator have voltage signals applied, the instantaneous responses are obtained as shown in Figure 2. The voltage signal is set as v input = 5sin(2πf ), with f being tested at 50 Hz, 80 Hz, and 200 Hz. The simulations of pinched hysteresis loops, referred to q M -v C characteristics, can be got from Multisim 12 software, which are shown in Figure 2.

Memcapacitor-Based Chaotic Oscillator
Based on the memcapacitor model in Equation (7), a chaotic oscillator is designed as shown in Figure 3, which contains conductances G1 and G2, capacitor C1, inductor L, and memcapacitor CM. where v1 = (a + b|σM(t)|)qM. Let the circuit parameters be chosen as shown in Table 1 with initial conditions (0.02, 0, 0, 0). Then, System (8) is chaotic with chaotic attractors as shown in Figure 4.

Memcapacitor-Based Chaotic Oscillator
Based on the memcapacitor model in Equation (7), a chaotic oscillator is designed as shown in Figure 3, which contains conductances G 1 and G 2 , capacitor C 1 , inductor L, and memcapacitor C M .

Memcapacitor-Based Chaotic Oscillator
Based on the memcapacitor model in Equation (7), a chaotic oscillator is designed as shown in Figure 3, which contains conductances G1 and G2, capacitor C1, inductor L, and memcapacitor CM. where v1 = (a + b|σM(t)|)qM. Let the circuit parameters be chosen as shown in Table 1 with initial conditions (0.02, 0, 0, 0). Then, System (8) is chaotic with chaotic attractors as shown in Figure 4. Taking the current i L through the inductor, the voltage v 2 across the capacitor and the charge q M on the memcapacitor as state variables, a set of four first-order state equations can be obtained, as follows: where v 1 = (a + b|σ M (t)|)q M . Let the circuit parameters be chosen as shown in Table 1 with initial conditions (0.02, 0, 0, 0). Then, System (8) is chaotic with chaotic attractors as shown in Figure 4.   Variable 0.5 The system Lyapunov exponents were LE1 = 0.0382, LE2 = 0.0067, LE3 = -0.0097, and LE4 = -0.2142. The Lyapunov dimension was DL = 3.1643. Figure 5a shows the Poincaré map on z = 0 and Figure 5b shows the Poincaré maps with initial condition v2(0) varying in the range of (-0.08, 0.08), showing that the system's dynamic behaviors were affected by initial conditions. All the above results verified that System (8) is chaotic [23].  Figure 5a shows the Poincaré map on z = 0 and Figure 5b shows the Poincaré maps with initial condition v 2 (0) varying in the range of (−0.08, 0.08), showing that the system's dynamic behaviors were affected by initial conditions. All the above results verified that System (8) is chaotic [23].

Equilibrium Points
The system equilibrium set can be calculated using E = {(x, y, z, w) | iL = v2 = qM = 0, σM = c} by solving the equations of = = = = 0, in which there is a real constant parameter c.
The Jacobian matrix J at this equilibrium set E is:

Equilibrium Points
The system equilibrium set can be calculated using The Jacobian matrix J at this equilibrium set E is: The corresponding characteristic equation is: It is obvious that the characteristic equation has one zero value and three nonzero values. Let Then, according to the Routh-Hurwitz condition, the system is stable if: When the parameters are set as in Table 1, one can find that: To make the equilibrium set E unstable, which is a necessary condition for the possible existence of chaos, the constant c should satisfy: |c| < 1.4002 or |c| > 2.7431 (13) Equation (13) demonstrates that the dynamical behaviors of the chaotic circuit (8) were strongly dependent on the memcapacitor internal state variable σ M . For example, if the system parameters were set as in Table 1, with c = 1, then four eigenvalues at the equilibrium set E were obtained as: (14) In this case, the equilibrium set E was unstable with one zero root, two complex conjugate roots with negative real parts, and one positive real root. Thus, a self-excited attractor could be generated via excitation from the unstable focal point in E.

Parameters Region
When the inductor L increased gradually with other circuit parameters fixed as in Table 1, the bifurcation diagram of the state variable i L is shown in Figure 6a, where the orbits of the system started from chaotic behavior and then entered into periodic behavior via the reverse period-doubling bifurcation route. After that, the orbits returned to chaotic behavior through the forward period-doubling bifurcation route, then jumped into chaotic behavior, and finally approached infinity in the range of L > 1. 25. The corresponding Lyapunov exponent spectra are presented in Figure 6b, where the maximum Lyapunov exponent was positive within chaotic regions and equalled zero within periodic regions. Several periodic windows can be observed in Figure 6b, which match well with that of the bifurcation diagram shown in Figure 6a. doubling bifurcation route. After that, the orbits returned to chaotic behavior through the forward period-doubling bifurcation route, then jumped into chaotic behavior, and finally approached infinity in the range of L > 1. 25. The corresponding Lyapunov exponent spectra are presented in Figure 6b, where the maximum Lyapunov exponent was positive within chaotic regions and equalled zero within periodic regions. Several periodic windows can be observed in Figure 6b, which match well with that of the bifurcation diagram shown in Figure 6a. The new system also had coexisting bifurcations. When the inductor G1 varied with other circuit parameters fixed as in Table 1, the bifurcation diagram of the state variable iL is shown in Figure 7a, where the red orbit started from initial conditions (0, -0.03, 0, 0), and the blue one starts from (0, 0.03, 0, 0). The symmetrically coexisting bifurcation orbits were generated by the symmetrical coexisting attractors, which are shown in Figure 9. The corresponding Lyapunov exponent spectrum is shown in Figure 7b. From the bifurcation diagram and the Lyapunov exponent spectrum, it can be seen that the system orbit started from a limit cycle and then turned into the chaotic state through perioddoubling bifurcations in the range of (0, 0.6). When the parameter G1 > 0.6, the bifurcation diagram and the Lyapunov exponent spectrum had blank areas since the system was diverging with no solution. Between the two blank areas, the system stayed in a periodic state.  The new system also had coexisting bifurcations. When the inductor G 1 varied with other circuit parameters fixed as in Table 1, the bifurcation diagram of the state variable i L is shown in Figure 7a, where the red orbit started from initial conditions (0, −0.03, 0, 0), and the blue one starts from (0, 0.03, 0, 0). The symmetrically coexisting bifurcation orbits were generated by the symmetrical coexisting attractors, which are shown in Figure 9. The corresponding Lyapunov exponent spectrum is shown in Figure 7b. From the bifurcation diagram and the Lyapunov exponent spectrum, it can be seen that the system orbit started from a limit cycle and then turned into the chaotic state through period-doubling bifurcations in the range of (0, 0.6). When the parameter G 1 > 0.6, the bifurcation diagram and the Lyapunov exponent spectrum had blank areas since the system was diverging with no solution. Between the two blank areas, the system stayed in a periodic state.
period-doubling bifurcation route, then jumped into chaotic behavior, and finally approached infinity in the range of L > 1. 25. The corresponding Lyapunov exponent spectra are presented in Figure 6b, where the maximum Lyapunov exponent was positive within chaotic regions and equalled zero within periodic regions. Several periodic windows can be observed in Figure 6b, which match well with that of the bifurcation diagram shown in Figure 6a. The new system also had coexisting bifurcations. When the inductor G1 varied with other circuit parameters fixed as in Table 1, the bifurcation diagram of the state variable iL is shown in Figure 7a, where the red orbit started from initial conditions (0, -0.03, 0, 0), and the blue one starts from (0, 0.03, 0, 0). The symmetrically coexisting bifurcation orbits were generated by the symmetrical coexisting attractors, which are shown in Figure 9. The corresponding Lyapunov exponent spectrum is shown in Figure 7b. From the bifurcation diagram and the Lyapunov exponent spectrum, it can be seen that the system orbit started from a limit cycle and then turned into the chaotic state through perioddoubling bifurcations in the range of (0, 0.6). When the parameter G1 > 0.6, the bifurcation diagram and the Lyapunov exponent spectrum had blank areas since the system was diverging with no solution. Between the two blank areas, the system stayed in a periodic state.

Similar Bifurcation Structures with Initial Conditions
Different bifurcation parameters usually lead to different bifurcation structures. However, there exist similar bifurcation structures with different initial conditions, which is rare compared with other chaotic systems. When we set the circuit parameters as given in Table 1, similar bifurcation diagrams with respect to i L (0), v 2 (0), and q M (0) were discovered, as presented in Figure 8. Although these diagrams depended on the various bifurcation parameters, it is remarkable that the bifurcation structures were almost the same and all symmetrical about the origin. As the bifurcation parameters increased gradually, the system orbit went to a chaotic status via period-doubling bifurcations. Then, the dynamics of the system settled down to periodic behaviors via reverse period-doubling bifurcations. The origin was the boundary of the two processes. The corresponding Lyapunov exponent spectra are given in Figure 8d-f, where the graphs also have similar shapes. diagrams depended on the various bifurcation parameters, it is remarkable that the bifurcation structures were almost the same and all symmetrical about the origin. As the bifurcation parameters increased gradually, the system orbit went to a chaotic status via period-doubling bifurcations. Then, the dynamics of the system settled down to periodic behaviors via reverse period-doubling bifurcations. The origin was the boundary of the two processes. The corresponding Lyapunov exponent spectra are given in Figure 8d-f, where the graphs also have similar shapes. Furthermore, the system's dynamic maps are used to illustrate the similar initial-conditiontriggered bifurcation structures, which are shown in Figure 9. These dynamic maps describe different dynamical regions with respect to the bifurcation parameters iL(0), v2(0), and qM(0), where the red areas indicate a chaotic field, the blue regions represent a periodic status, and yellow areas show unbounded zones. It is obvious the two dynamic maps own a similar distribution structure. Thus, we can conclude that the initial conditions iL(0), v2(0), and qM(0) had a similar dynamic influence for the presented system in the especial parameter spaces, which is not common in the other chaotic systems. Furthermore, the system's dynamic maps are used to illustrate the similar initial-conditiontriggered bifurcation structures, which are shown in Figure 9. These dynamic maps describe different dynamical regions with respect to the bifurcation parameters i L (0), v 2 (0), and q M (0), where the red areas indicate a chaotic field, the blue regions represent a periodic status, and yellow areas show unbounded zones. It is obvious the two dynamic maps own a similar distribution structure. Thus, we can conclude that the initial conditions i L (0), v 2 (0), and q M (0) had a similar dynamic influence for the presented system in the especial parameter spaces, which is not common in the other chaotic systems.

Extreme Multistability and Coexisting Attractors
Multistability is a common phenomenon in many nonlinear dynamical systems, corresponding to the coexistence of more than one stable attractor for the same set of system parameters [24]. When infinitely many attractors coexist for the same set of system parameters, multistability is referred to as extreme multistability [25]. The previously published literature have reported that the dynamical stability of memristive systems are heavily dependent on the initial conditions, which easily leads the system to generate multistability or even extreme multistability [26][27][28][29].
One of the main features of extreme multistability is that the system track can present bifurcation without varying any system parameter. When the circuit parameters were set as in Table 1, with initial condition iL(0) varying, the resulting bifurcation diagram of the state variable iL is shown in Figure 10a, where the system presents extreme multistability. It is remarkable to see that the bifurcation diagram is symmetrical about the origin, which came from symmetrical coexisting

Extreme Multistability and Coexisting Attractors
Multistability is a common phenomenon in many nonlinear dynamical systems, corresponding to the coexistence of more than one stable attractor for the same set of system parameters [24]. When infinitely many attractors coexist for the same set of system parameters, multistability is referred to as extreme multistability [25]. The previously published literature have reported that the dynamical stability of memristive systems are heavily dependent on the initial conditions, which easily leads the system to generate multistability or even extreme multistability [26][27][28][29].
One of the main features of extreme multistability is that the system track can present bifurcation without varying any system parameter. When the circuit parameters were set as in Table 1, with initial condition i L (0) varying, the resulting bifurcation diagram of the state variable i L is shown in Figure 10a, where the system presents extreme multistability. It is remarkable to see that the bifurcation diagram is symmetrical about the origin, which came from symmetrical coexisting attractors. As initial condition i L (0) increased gradually within the region of [−0.94, 0], the system orbit started from a limit cycle and turned into a chaotic state through period-doubling bifurcations, showing several periodic windows. Within the region of [0, 0.94], the system orbit settled down to periodic behavior from the chaotic state through reverse period-doubling bifurcations, which is the reverse evolutionary process of that in the region [−0.94, 0]. The corresponding Lyapunov exponent spectra are shown in Figure 10b, where the maximum Lyapunov exponents stay zero with limit cycles but were positive in chaotic states, consistent with the bifurcation diagram.

Extreme Multistability and Coexisting Attractors
Multistability is a common phenomenon in many nonlinear dynamical systems, corresponding to the coexistence of more than one stable attractor for the same set of system parameters [24]. When infinitely many attractors coexist for the same set of system parameters, multistability is referred to as extreme multistability [25]. The previously published literature have reported that the dynamical stability of memristive systems are heavily dependent on the initial conditions, which easily leads the system to generate multistability or even extreme multistability [26][27][28][29].
One of the main features of extreme multistability is that the system track can present bifurcation without varying any system parameter. When the circuit parameters were set as in Table 1, with initial condition iL(0) varying, the resulting bifurcation diagram of the state variable iL is shown in Figure 10a, where the system presents extreme multistability. It is remarkable to see that the bifurcation diagram is symmetrical about the origin, which came from symmetrical coexisting attractors. As initial condition iL(0) increased gradually within the region of [−0.94, 0], the system orbit started from a limit cycle and turned into a chaotic state through period-doubling bifurcations, showing several periodic windows. Within the region of [0, 0.94], the system orbit settled down to periodic behavior from the chaotic state through reverse period-doubling bifurcations, which is the reverse evolutionary process of that in the region [−0.94, 0]. The corresponding Lyapunov exponent spectra are shown in Figure 10b, where the maximum Lyapunov exponents stay zero with limit cycles but were positive in chaotic states, consistent with the bifurcation diagram. To show more details, various typical coexisting attractors are displayed in Figure 11. When the circuit parameters were set as in Table 1, with different initial conditions, the main coexisting regimes were symmetric pairs of limit cycles, chaotic attractors, and point attractors. Figure 11a  To show more details, various typical coexisting attractors are displayed in Figure 11. When the circuit parameters were set as in Table 1, with different initial conditions, the main coexisting regimes were symmetric pairs of limit cycles, chaotic attractors, and point attractors. Figure 11a shows limit cycles coexisting with attractors for initial conditions (0.02, ±0.06, 0, 0). Figure 11b,c shows coexisting chaotic attractors with initial values (0, ±0.04, 0, 0) and (0.02, 0, 0, 0), respectively. Figure 11d displays a point attractor with initial conditions (−1, 0.18, 0, 0). Since the system track presented bifurcations with initial conditions and had infinitely many coexisting attractors for the same set of system parameters, the system presented typical extreme multistability.
Basins of attraction can clearly display the distribution of different coexisting attractors. As shown in Figure 11, the system had six kinds of coexisting attractors. When initial conditions i L (0) and v 2 (0) were set as variables, with q M (0) = 0 and σ M (0) = 0, the corresponding basin of attraction is displayed in Figure 12a, where coexisting attractor distributions are painted with different colors. As shown in Figure 12a, the system was divergent with no attractor in most blank areas and generated a narrow distribution of each kind of the coexisting attractors in the middle areas. The evolution process of coexisting attractors is shown in Figure 12b, which contains a complete process of period-doubling bifurcations and reverse period-doubling bifurcations as v 2 (0) decreases. The system orbit started from a point attractor (Type 6) to a limit cycle (Type 1), and then turned to a chaotic attractor (Type 5) via period-doubling bifurcations (Type 3). As v 2 (0) decreased further, the chaotic attractor (Type 5) returned to a limit cycle (Type 2) via reverse period-doubling bifurcations (Type 4) and finally settled into a point attractor (Type 6).  Figure 11d displays a point attractor with initial conditions (-1, 0.18, 0, 0). Since the system track presented bifurcations with initial conditions and had infinitely many coexisting attractors for the same set of system parameters, the system presented typical extreme multistability. Basins of attraction can clearly display the distribution of different coexisting attractors. As shown in Figure 11, the system had six kinds of coexisting attractors. When initial conditions iL(0) and v2(0) were set as variables, with qM(0) = 0 and σM(0) = 0, the corresponding basin of attraction is displayed in Figure 12a, where coexisting attractor distributions are painted with different colors. As shown in Figure 12a, the system was divergent with no attractor in most blank areas and generated a narrow distribution of each kind of the coexisting attractors in the middle areas. The evolution process of coexisting attractors is shown in Figure 12b, which contains a complete process of perioddoubling bifurcations and reverse period-doubling bifurcations as v2(0) decreases. The system orbit started from a point attractor (Type 6) to a limit cycle (Type 1), and then turned to a chaotic attractor (Type 5) via period-doubling bifurcations (Type 3). As v2(0) decreased further, the chaotic attractor (Type 5) returned to a limit cycle (Type 2) via reverse period-doubling bifurcations (Type 4) and finally settled into a point attractor (Type 6).
Similarly, the basin of attraction with respect to initial conditions qM(0) and σM(0) is displayed in Figure 12c, which is roughly symmetric about the origin and each attractor distribution appears in a narrow reverse 'S' shape. The corresponding evolution process of coexisting attractors is shown in Figure 12d. Similarly, the basin of attraction with respect to initial conditions q M (0) and σ M (0) is displayed in Figure 12c, which is roughly symmetric about the origin and each attractor distribution appears in a narrow reverse 'S' shape. The corresponding evolution process of coexisting attractors is shown in Figure 12d.

Experimental Results
An analog electronic circuit was built to physically realize the above-presented chaotic oscillator to verify the basic dynamic behaviors of the new system.
Since the inductor and the capacitor in the system were not standard, which made the circuit difficult to design and implement, an equivalent circuit was designed to realize the system. The circuit schematic of the experimental circuit is shown in Figure 13, and the corresponding chaotic attractors obtained by Multisim simulation are displayed in Figure 14.
In practical experiments, since there are tolerances in resistors and capacitors, it is necessary to adjust the actual resistance values in the analog circuit, which will lead to some deviations between experimental results and the simulation ones. The experimental results shown on the oscilloscope are depicted in Figure 15. The chip AD633JN was chosen as the analog multiplier and LF347N as the operational amplifier with reference voltages of ±15 V. It is clear that the dynamical behaviors observed from the experimental circuit were generally similar with those displayed via numerical simulations.

Experimental Results
An analog electronic circuit was built to physically realize the above-presented chaotic oscillator to verify the basic dynamic behaviors of the new system.
Since the inductor and the capacitor in the system were not standard, which made the circuit difficult to design and implement, an equivalent circuit was designed to realize the system. The circuit schematic of the experimental circuit is shown in Figure 13, and the corresponding chaotic attractors obtained by Multisim simulation are displayed in Figure 14.   In practical experiments, since there are tolerances in resistors and capacitors, it is necessary to adjust the actual resistance values in the analog circuit, which will lead to some deviations between experimental results and the simulation ones. The experimental results shown on the oscilloscope are depicted in Figure 15. The chip AD633JN was chosen as the analog multiplier and LF347N as the operational amplifier with reference voltages of ±15 V. It is clear that the dynamical behaviors observed from the experimental circuit were generally similar with those displayed via numerical simulations.

Conclusions
In this paper, a new memcapacitor model and a novel memcapacitor-based chaotic circuit are presented. The system extreme multistability was analyzed, including bifurcation diagrams, Lyapunov spectra, coexisting attractors, coexisting bifurcations, and basins of attraction of various

Conclusions
In this paper, a new memcapacitor model and a novel memcapacitor-based chaotic circuit are presented. The system extreme multistability was analyzed, including bifurcation diagrams, Lyapunov spectra, coexisting attractors, coexisting bifurcations, and basins of attraction of various attractors. Moreover, the new memcapacitor-based system was realized using an experimental circuit, which agreed well with the numerical simulations and verified the theoretical analysis results. Due to the rich and unusual complex dynamical characteristics of the proposed memcapacitor system, it was deemed that it would find some novel and non-traditional applications in engineering and technology in the future. In our future works, we will continue to try to build physical memcapacitors and explore special dynamics in memcapacitive circuits.