All Pinched Hysteresis Loops Generated by (α, β) Elements: in What Coordinates They May be Observable

Two existing theorems for studying pinched hysteresis loops generated by nonlinear higher-order elements from Chua’s table are reformulated, namely the generalized homothety theorem and the associated Loop Location Rule, specifying the coordinates where the hysteresis may occur, and the <inline-formula> <tex-math notation="LaTeX">$\omega ^{2}$ </tex-math></inline-formula> criterion theorem for computing the corresponding loop areas. It is demonstrated in this work that the pinched hysteresis loops are also generated in other coordinates than those predicted by the Loop Location Rule, and all these possible coordinates are found. The <inline-formula> <tex-math notation="LaTeX">$\omega ^{2}$ </tex-math></inline-formula> criterion is generalized to computing the areas of all hysteresis loops that may be observable.


I. INTRODUCTION
Nonlinear (α, β) two-terminal elements, also referred to as the higher-order elements, organized in Chua's table of fundamental electrical elements [1], are an important tool for the so-called predictive modeling of complex nonlinear dynamical systems and processes including the phenomena in molecular and nanoscale devices [2]. The (α, β) element preserves an unambiguous constitutive relation between the pair of quantities v (α) and i (β) , called constitutive variables, where the positive, zero, or negative integers α and β are orders of the time derivatives or integrals of terminal voltages and currents of the element. The most widely known electrical elements from the infinite set of (α, β) elements are the resistors, capacitors, inductors, and memristors with the coordinates (0,0), (0,-1), (-1,0), and (-1,-1). The synthesis of the model of a concrete system is based on the selection of the (α, β) elements with the relevant nonlinear constitutive relations and on their proper interconnection. This approach facilitates the physical insight into the mechanism of complex phenomena (see for example the predictive model of the The associate editor coordinating the review of this manuscript and approving it for publication was Cihun-Siyong Gong . Josephson junction in [2]) and also significantly decreases the simulation times (see for example the simulation of predictive models of the spin-torque nano-oscillator in [3]).
The substantial feature a model built in such a way is its above-mentioned predictability, which means the ability to predict the behavior of the modeled object under general conditions of its interaction with the surroundings, thus, among others, for arbitrary types of driving signals, etc. An example of the predictive model of a nonlinear capacitor is its coulomb -volt characteristic, whereas its small-signal impedance as a function of frequency is a typical non-predictive model. The latter model, which is effective for a concrete operating point, cannot be used for predicting the behavior of the capacitor under its general (large-signal) interaction with the remainder of the circuit.
The predictability of a model made up of nonlinear (α, β) elements starts from the thesis that the predictable model must consist of building blocks whose models are predictable in their own right.
It turns out that the concept of the (α, β) two-terminal elements can also be used for effective modeling of systems outside of electrical engineering, where, instead of the voltage and current and their derivatives and integrals, VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the pair of effort (e) and flow ( ) [4] variables, typical for a given technical area, is used. Let us mention force and velocity for mechanical engineering, temperature and entropy flow for thermodynamics, chemical potential and molar flow for chemical engineering, etc. The works [5]- [7] deal with special mechanical elements, which are, however, the mechanical equivalents of previously introduced electrical memristors and memcapacitors [8], [9]. Also the so-called inerter, a newly discovered mechanical element [10], has its existing archetype in Chua's classical table of (α, β) electrical elements [11]. The procedures of predictive modeling, developed in electrical engineering, can be therefore applied in a multi-disciplinary way with the utilization of the relevant version of the table of (α, β) elements. It must be obtained via transforming Chua's table from electrical engineering into a given scientific branch via respective analogies, for example the well-known electric-mechanical voltage-force or voltagevelocity analogy.
The above modeling approach can be advantageous for studying interesting hysteretic effects, which have been experimentally observed and intensively studied in various branches of science. The hysteresis loops, pinched at the origin of the coordinates of appropriate measured variables, will be denoted hereinafter as the PHL (pinched hysteresis loops). Their shape and area depend on the parameters of the system and the driving signal. There have been earlier works dealing with hysteresis in the current-voltage characteristics of the electric arc or incandescent lamps [12]. The hysteretic effects in the i-v characteristics of the CdSe point contact diodes [13] were published in the same year the memristor was introduced into the circuit theory as a hypothetical element [14], and 37 years before the manufacture of a nanodevice that exhibited the features of the memristor [15] as the (α, β) element for α = β = −1. The pinched hysteresis that accompanies the deformations of inelastic materials and the cyclic stress of buildings simulating an earthquake are documented in [16]. Works from the area of organic and inorganic nature model the memory effects in gases [17], in cylindrical protein polymers (microtubules) [18], in nanofluidic [19] and quantum [20] memristors, in synaptic junctions [21], in hafnium oxide-based ferroelectrics [22], in solutions with measuring electrodes [23], in muscle fibers [24], in human skin [25], or in the behavior of plants and primitive organisms [26], [27] with the utilization of the (-1,-1) element, i.e. the memristor. Novel pinched nonlinear hysteretic structural models for effective processing of measured hysteresis loops in various types of complex systems are suggested in [28].
In the past, the causes of the studied hysteretic behavior were often identified incorrectly. The respective systems are denoted in [2] by the attribute ''mistaken identity''. A typical example is the original Huxley-Hodgkin axon model [29], containing time-varying conductances whose manifestations were interpreted, in agreement with [30], as inductive reactances. Later on, this non-predictive model was substituted by correct models of the Potassium and Sodium ion channels, which were identified as memristors [31].
The clue for revealing the true identity of systems that generate the PHLs can be the Generalized Homothety Theorem (GHT). In [32], it is defined as follows: The movement of the operating point in the (v (α) , i (β) ) space along the nonlinear constitutive relation of an (α, β) element from Chua's table is accompanied by the movement of the operating point in the (v (α+1) , i (β+1) ) space along a pinched hysteresis loop. The loops that correspond to the driving signals of the same levels but various frequencies are homothetic entities with the homothetic center at the origin of the (v (α+1) , i (β+1) ) space, and the scale factor equal to the ratio of frequencies. The areas within the lobes of the loops increase with the square of the frequency.
The GHT, applicable to (α, β) elements with arbitrary waveforms v (α) and i (β) , specifies three regularities, associated with this element: The LLR can have useful consequences. For example, if the PHLs are observed in voltage-current, thus v (0) -i (0) coordinates, then the source of such a behavior can be the (-1,-1) element, which is the memristor defined by the constitutive relation between the v (−1) -i (−1) variables (time integrals of the voltage and current, denoted as the flux and charge). For the PHLs studied in the stress-strain characteristics of cyclically stressed polymers, the (-1,-2) mechanical element, whose constitutive variables are the integral of force vs. the second integral of velocity, can be responsible for the hysteretic behavior: Considering that the effort is the force and the flow is the velocity, the stress-strain coordinates can be interpreted as a force-displacement pair (the first integral of the velocity). The corresponding electrical (-1,-2) element is the memcapacitor.
To verify the hypothesis offered by the LLR, namely that a concrete (α, β) element is responsible for the hysteretic behavior, the HR, FR, and other auxiliary rules dealing with the symmetry and frequency dependence of the loop area can be used.
According to [33], for a sinusoidal driving signal, the loop is formed by two lobes centrally symmetric with respect to the origin of the coordinates. Many papers have been devoted to evaluating the area of the lobes of PHLs generated by various mem-systems, and how these areas depend on the frequency [34]- [43]. According to [34], the lobe area S can be found from the area Å between the curve of the constitutive relation g() and the chord linking two terminal operating points (see Fig. 1) by using the formula S = ω 2 Å, which will be denoted below as the ω 2 criterion. This criterion is derived in [34] for a memristor driven by a charge with harmonic waveform. The derivation is made from more general regularities revealed in [34], namely that the area within the loop is directly proportional to the value of an exactly defined part of the (co)action potential. The ω 2 criterion is closely related to the FR, which is a part of the homothety theorem. FR exactly describes the evolution of the loop area if the frequency of a general driving signal is changed, but it does not provide a formula for computing this area. On the other hand, the ω 2 criterion enables evaluating the area, but only for harmonic driving signals.
The cause of the asymmetry of the PHL can be seen in that the corresponding system cannot be properly modeled only by one ideal (α, β) element but via a more complex scheme of interconnected models, which in their entirety exhibit a more complicated behavior [27], [44]. As one more reason of the asymmetry, the model of some systems cannot be compounded from a mere set of two-terminal (α, β) elements [3]. A recommended guideline for building such models is the analysis of physical processes in the system with the aim to obtain their equations of motion. These equations are the starting point for the synthesis of the predictive model via procedures described in [3], [45]. Also the above information about the frequency dependence of the area of generated PHL can be useful for building up the model. An instructive example is given in the work [17], dealing with the analysis of the hysteretic relationship between the discharge current and the gap voltage in helium dielectric barrier discharges.

II. MOTIVATION, CONTRIBUTIONS, AND PAPER ORGANIZATION
Based on the state-of-the-art analysis in Section I, the following two questions come into consideration: 1. Is there also some other pair of variables in addition to the first derivatives of the constitutive variables for which, according to the LLR, the PHLs are observable in their coordinates? 2.
If yes, what are the corresponding regularities? Can the ω 2 criterion also be generalized to computing the areas of these loops?
As will be shown in Section III, the answer to the question 1 is yes. As a consequence, the LLR cannot be directly used for an unambiguous identification of the element which is responsible for the hysteretic behavior. If the hysteresis can be observable in various coordinates, it is not obvious which of them can be used for specifying the coordinates of the element in the table. The other consequence is the necessity of engaging in the problems from question 2.
The motivation for this work is to resolve all the above issues. The remainder of the paper is organized as follows. The LLR is reformulated in Section III such that all possible coordinates are found where the PHLs generated by an arbitrary (α, β) element can be observable. Section IV deals with the computation of PHL areas. The existing ω 2 criterion is extended, and the relation between the areas denoted as + and -in Fig. 1 and this criterion is clarified. Simulations in Section V demonstrate the usefulness of the above new pieces of knowledge in electrical and mechanical engineering. Finally, Section VI summarizes the results and suggests possible ways of continuing in the research.
To generalize the notation, the effort (e) and flow ( ) variables will be considered hereinafter instead of electrical voltage (v) and current (i). In Section V, describing the simulations of mechanical and electrical systems, the effort and flow variables will be specified for concrete domains.

III. MODEL OF COMPLETE SYSTEM OF PHLS IN CHUA'S TABLE: REVISITING THE LLR
Generally, the PHL appears in the e (α+j) − f (β+k) plane, j, k = 0, if and only if both variables simultaneously cross the zero value. Let the (α, β) element be driven by a biased waveform where u 0 provides a proper operating point of the element, and U , ω are the amplitude and angular frequency. The oddand even-order derivatives of the signal (1) are Faà di Bruno's formula [46] for the n th -order derivative of the quantity f (β) = g(e (α) ) holds: The symbol g () denotes a higher-order derivative of the constitutive relation with respect to e (α) . The summation is performed via all n-tuple nonnegative integers (m 1 , ..m n ), VOLUME 8, 2020 which are the solution to the Diophantine equation Even-order derivatives of f (β) cannot form the PHL with either odd-or even-order derivatives of e (α) . To generate a loop, the right side of (3) for n even would have to contain either sine or cosine components in each of its addends. It would enable the generation of PHL in combination with the odd or even order of derivative (2). However, this is not fulfilled, since two of the possible solutions to (4)  For an odd-order derivative of f (β) , the PHL cannot appear for the even-order derivative of e (α) , since with n being odd, all solutions to the Diophantine equation (4) contain nonzero parts m j with odd indices j. As a result, the sine component contained in every addend (3) is decisive for crossing odd-order time derivative of f (β) by zero value.
When driving the (α, β) element via a harmonic signal, the PHL is, for reasons given above, generated only in the e (α+j) − f (β+k) planes, where j and k are odd positive integers. The layout is illustrated in Fig. 2. The symbols of the loop or cross placed in the coordinates (α + j, β + k), j, k = 1, 2, 3, .. inform about whether the PHL is or is not generated in the e (α+j) -f (β+k) planes.

IV. LOOP AREA COMPUTATION: REVISITING THE ω 2 CRITERION
The area of the PHL lobe can be computed as a generalized content C or co-content C * [47] generated within one halfperiod of the sine component, when the complete PHL lobe was formed. For a general pair of signals in the e (α+j) -f (β+k) plane, the content and co-content are Integrating by parts the definition equations (5) leads to regularities which hold for every j, k: The last relation in (6) is the Legendre transformation between the content and co-content. With j, k being odd, the PHLs are generated, and the content and co-content differ only in the sign, because e (α+j) is equal to zero at the beginning and end of the first half-period.
Eq. (1) and (2) yield the algebraic relation between e (α) and e (α+2) , namely e (α+2) . The area of the loop lobe in the e (α+1) -f (β+1) plane will be The integral in (7) is the content of the constitutive relatioñ g(), which is transformed from the original relation g() by moving the origin of the coordinates into the operating point (u 0 , y 0 ). This content corresponds to the negative algebraic sum of the areas denoted in Fig. 1 as + and -.
Applying the rule (6) yields a result which is analogous to (7) Eq. (8) is the classical ω 2 criterion [34] applied to common (α, β) elements. Relations (6) and (2) lead to the following ω 2 criterion, which holds for arbitrary PHLs: It follows from (2) that e (α+j+2) = -ω 2 e (α+j) for every j ≥ 0. Considering the definitions (5) and the fact that the content and co-content differ only in signs for odd j, k, it must hold for the effort-type excitation odd j, k : C * j+2,k = ω 2 C j,k Eq. (6) and (9) can be used for computing the area of every PHL represented by the content C j,j for an arbitrary odd j. An example is the area of the loop generated in the e (α+3) - where G is a difference in the generalized conductances of the element between the end and beginning of the half-period of excitation. Geometrically, it is the difference between the tangents of the angle at which the PHL is leaving the f (β+1)e (α+1) origin and the angle at which the PHL is coming back.
The value of C j,j for an arbitrary j > 0 is the starting point for determining the areas of the PHLs whose symbols in Fig. 2 appear on the horizontal line which corresponds to the value j. The areas are determined according to the relation (10).

V. SIMULATIONS A. MECHANICAL INERTER
The vibration absorber in Fig. 3 containing a linear inerter is designed in [10]. It consists of a pair of springs with stiffness k 1 and k 2 , a damper with the damping coefficient c, and an inerter with the inertance b, which for the linear inerter is the constant of proportionality between the acceleration a and the applied force F. Then the corresponding constitutive relation of the linear inerter is where v is the velocity. Note that the inerter is therefore the (0,1) element in the framework (effort, flow) = (force, velocity). The absorber is designed as a mechanical bandpass rejection filter such that the 10 Hz vibrations, causing the deflection z, are not projected onto the coordinate x of the body with the mass M. The parameters of the absorber are as follows [10]: M = 10 kg, k 1 = 9 × 10 4 Nm −1 , The results of a transient analysis of the absorber driven into the node z by a 10 Hz sinusoidal vibration with a velocity amplitude of 0.1 ms −1 are available in [11]. It is shown therein that the prospective nonlinearity, introduced into the constitutive relation (12), makes the filtering of the vibration worse. The nonlinearity considered in [11] consists in the limitation of the acceleration to the levels ±100/b if the absolute value of the force exceeds the value of 100 N. However, this means that the first derivative of the acceleration is not a continuous function of time. In order to utilize this mechanical system for verifying the Loop Location Rule, i.e. that the PHLs appear in the coordinates of the appropriate odd derivatives of constitutive variables, this limitation must be modeled via a function with continuous derivatives, for example where the parameter m determines the limit value of the acceleration. For the inerter driven by a small force F →0, the nonlinear constitutive relation (13) passes into the linear form (12). For a large excitation force, the acceleration is limited to the value m/b, the same as in [11]. The parameter m = 50 is chosen for the subsequent simulations with a view to accentuating the nonlinear effects. Fig. 3 (b) demonstrates the absorber transition to the steady state. It is shown that the nonlinearity of the inerter is responsible for the acceleration limits m/b = 2.193 ms −2 . According to the conclusions from Section III, the PHLs of the sine-driven inerter should be drawn in coordinates of the first derivatives of the force and acceleration. However, the inerter is not driven by a harmonic signal during the transient in Fig. 3 (b). Consequently, a set of PHLs may be observable in Fig. 3 (b), which approaches to a steady-state hysteretic pattern.
The results of the steady-state analysis are summarized in Fig. 4. Since the computation of higher-order derivatives from the simulated waveforms is subject to nonacceptable errors, the patterns in Fig. 4 were obtained from analytic formulae of these derivatives, derived from (13) and on the assumption of sinusoidal force acting on the inerter with the amplitude extracted from the original steady-state simulation. Fig. 4 confirms the LLR, since the PHLs appear in coordinates that are also displayed in Fig. 2. The symmetries of some characteristics, which can be traced up from the simulation results, are particularly remarkable.

B. NONLINEAR RESISTIVE CIRCUIT WITH SILICON DIODES
Consider the electrical resistive (0,0) element as two identical antiparallel diodes with the resulting hyperbolic-sine-type current-voltage constitutive relation where I s , n, and v T are the saturation current, emission coefficient, and thermal voltage, respectively. In the simula-VOLUME 8, 2020 where R 0 , R off , and R on are the initial, OFF, and ON resistances of the memristor, Q 0 is a charge necessary for switching between the ON and OFF states. The memristor is driven by the signal (1), namely by the charge q = i (−1) with the  amplitude U ≡ Q max and the offset u 0 ≡ q 0 . Considering (2) with k = 1, the excitation can be accomplished via sinusoidal current with the frequency-dependent amplitude I max = ωQ max . The offset q 0 , specifying the initial operating point on the constitutive relation ϕ = f (q), is given by the initial value of the memristance R 0 . The loop graphs are again properly scaled. The maximum derivatives of the fluxes (k) max for k = 1 and 3 are 877.9 mV and 60.375 Vs −2 .
It follows from the duality principle for the (α, β) elements [49] that the ω 2 criterion (9) for the voltage excitation can be reworded for the current excitation via a mere interchange of the contents and co-contents. For a current-driven memristor with a general nonlinear constitutive relation, the area of the loop in the ϕ (3)q (3) plane will be given by the content Arranging (16) yields An interesting interpretation of (17) is presented in Fig. 7. The area of the triangle 0AB formed by the perpendicular crossing through the point I max and by two tangents to the PHL lobe, led by the vi origin, is It is obvious that, if the constitutive relation is concave, the lobe is clockwise oriented, thus C 1,1 > 0 and R < 0; for the convex constitutive relation, the opposite is true. In any case, the parentheses in (17) contain the sum of two terms with the same signs. The first term represents the area of the vi lobe, the second term is twice the area of the triangle circumscribed around this lobe.
It is proved in [34] for the HP memristor that the area of the vi lobe is equal to two-thirds of the S 0AB area. Connecting this piece of knowledge with (17) reveals the following regularity concerning the area of the v (2)i (2) lobe:

VI. CONCUSION
New pieces of knowledge presented in this work can be summarized as follows: 1. An arbitrary (α, β) element driven by a harmonic signal generates pinched hysteresis loops in the coordinates of odd-order time derivatives of its constitutive variables. The cause of the hysteresis is the nonlinearity of the constitutive relation of the element. 2. The algorithm of computing the area of the lobe of an arbitrary pinched hysteresis loop from item 1 is found. The relations between the areas of various pinched hysteresis loops are governed by the generalized ω 2 criterion, which allows an area computation based directly on the constitutive relation of the element. This criterion introduces order into researching the frequency dependence of higher-order hysteresis in nonlinear circuits. 3. The regularities 1 and 2 hold for arbitrary (α, β) elements, thus not only for (mem)resistors, (mem) capacitors and (mem)inductors. It is worth noting that they take effect independently of the type of the nonlinearity of the constitutive relation of the element. 4. The findings 1 and 2 hold for the (α, β) elements with a general pair of (effort, flow) constitutive variables. It makes them useful not only in the classical electrical engineering but also in other branches of science.
The above findings open up new possibilities of ongoing research. The following open issues are of particular interest: 1. Is it possible to generalize the Homothety Rule (HR, the part of The Generalized Homothety Theorem) to all loops according to Fig. 2? 2. What regularities and relationships govern the various forms of the characteristics (e (α+k) , f (β+l) ), k, l ∈ N+ (see the example of the characteristics of the inerter in Fig. 4)? 3. How does the loop area depend on frequency for a fixed amplitude of the driving variable e (α+k) , k ∈ N+? Hitherto known are only the results for the loop corresponding to the content C 1,1 for the memristor, k = 0 (the area increases with the square of frequency) and k = 1 (well-known frequency criterion: the area disappears with increasing frequency). 4. Searching for the identity of the element via the PHL analysis with the help of item 3. An exemplary work is [17] about researching the discharges in helium. The area of the PHL lobe grows with frequency, and the element identity is not currently known. 5. Researching into a complete system of PHLs (Loop Location Rule) for general waveforms of the constitutive variables of the element.
ZDENĚK BIOLEK received the Ph.D. degree in electronics and informatics from the Brno University of Technology, Czech Republic, in 2001. He is currently with the Department of Microelectronics, Brno University of Technology, and with the Department of EE, University of Defence, Brno, Czech Republic. Until the year 1993, he worked as an Independent Researcher with semiconductor company TESLA Rožnov. He is the author of unique electronic instruments associated with IC production and testing. He is also the author of several papers from the area of the utilization of variational principles in electrical engineering, and also from the field of memristors and mem-systems. He is the co-author of two books about memristive systems and modeling and simulation of special electronic circuits including switched-capacitor filters, switched dc-dc converters, and memristors. He is currently with the Department of EE, University of Defence (UDB), Brno, and the Department of Microelectronics, BUT. His scientific activity is directed to the areas of general circuit theory, frequency filters, mem-systems, and computer simulation of electronic systems. He is currently a Professor with BUT and UDB in the field of theoretical electrical engineering. He is a member of the CAS/COM Czech National Group of the IEEE. He has been a member of editorial boards of international journals including the International Journal of Electronics and Communications (AEU) and Electronics Letters.
VIERA BIOLKOVÁ (Member, IEEE) received the M.Sc. degree in electrical engineering from the Brno University of Technology, Czech Republic, in 1983. She joined the Department of Radio Electronics, in 1985, where she is currently working as a Research Assistant with the Department of Radio Electronics, Brno University of Technology, Czech Republic. Her research and educational interests include modeling of large-scale systems, signal theory, analog signal processing, memristors and memristive systems, optoelectronics, and digital electronics.
ZDENĚK KOLKA (Member, IEEE) received the M.Sc. and Ph.D. degrees in electrical engineering from the Brno University of Technology (BUT), Czech Republic, in 1992 and 1997, respectively. In 1995, he joined the Department of Radio Electronics, Brno University of Technology. He is currently a Professor of radio electronics with BUT. His scientific activity is directed to the areas of general circuit theory, computer simulation of electronic systems, and digital circuits. For years, he has been engaged in algorithms of the symbolic and numerical computer analysis of electronic circuits. He has published over 100 articles.