Stability Analysis of the Synchronous Solution in Arrays of Memristive Chua’s Circuits

This brief studies synchronization in an array of diffusively-coupled identical memristor Chua’s circuits (MCCs) with Dirichlet’s boundary conditions. We start from the description of MCCs in the flux-charge domain and the state space decomposition in invariant manifolds and then we apply three different techniques, of numerical and analytic nature, for finding the synchronization condition in terms of a circuit parameter. Two of these techniques provide conditions to guarantee synchronization on any manifold and for any attractor on a manifold, while a third technique allows analyzing synchronization for each specific attractor on a fixed manifold through more precise conditions. The pros and cons of the three techniques are discussed and an application to a more realistic memristor model is provided.

Abstract-This brief studies synchronization in an array of diffusively-coupled identical memristor Chua's circuits (MCCs) with Dirichlet's boundary conditions. We start from the description of MCCs in the flux-charge domain and the state space decomposition in invariant manifolds and then we apply three different techniques, of numerical and analytic nature, for finding the synchronization condition in terms of a circuit parameter. Two of these techniques provide conditions to guarantee synchronization on any manifold and for any attractor on a manifold, while a third technique allows analyzing synchronization for each specific attractor on a fixed manifold through more precise conditions. The pros and cons of the three techniques are discussed and an application to a more realistic memristor model is provided.

I. INTRODUCTION
N ETWORKS of coupled systems are receiving increasing attention from the scientific community, from time to time focusing on the topology, the node dynamics, the connections, the synchronization (global or partial), or other aspects [1], [2], [3], [4], [5], [6], [7]. In [8], arrays of diffusively coupled identical memristor Chua's circuits (MCCs) with Dirichlet's boundary conditions were analyzed. By using the flux-charge analysis method [9], [10], [11], it was shown that the state space can be decomposed in infinitely many invariant manifolds, where each manifold is characterized by an index depending upon the initial conditions in the voltagecurrent domain. On the manifold with index 0, the dynamics are analogous to that of a memristor-less array of Chua's circuits, while on a manifold with an index different from 0, we obtain the dynamics of an array of forced Chua's circuits with a forcing term depending on the index. Therefore, MCCs display an extremely rich dynamical scenario where, by changing the index (i.e., the initial conditions), different types of attractors (equilibrium points, limit cycles, chaotic attractors, etc.) Manuscript  and even multiple attractors for a given manifold [12] can coexist.
The analysis carried out in [8] is a first step in the understanding of complex phenomena displayed by MCCs. Yet, in that paper synchronization is analyzed only on the manifold with index 0 and only by numerically simulating the state equations. The main goal of this brief is to study the synchronization of MCCs also on manifolds with an index different from 0, by using techniques of both analytical and numerical nature.
The considered MCCs array is Laplacian and with diffusive connections, which ensures the existence of a minimal conductance value of the connecting resistors, beyond which the array synchronization is guaranteed. In this brief, we study synchronization conditions for an array of N diffusively-coupled identical MCCs through three different techniques, each one with its own pros and cons. The first one is the analytical technique proposed in [13], which provides a (sufficient) condition ensuring synchronization for any attractor on any manifold. The second one is a numerical method, based on the master stability function approach [14], which provides the synchronization condition for each attractor on each manifold. The third technique is a mixed analytic/numerical method, which generalizes the technique proposed in [13] to provide a (sufficient) condition ensuring synchronization for each specific attractor on a given manifold.
The results, obtained for N = 4, show that all conditions are easily scalable to arrays with a generic number N of MCCs and allow us to discuss the strengths and weaknesses of each technique. Moreover, we show via an example that the methods can be applied also to circuits containing more realistic memristor models.

II. ARRAY OF MEMRISTOR CHUA'S CIRCUITS
We analyze the global synchronization in a 1D array of N diffusively-coupled identical memristor Chua's circuits (MCCs). Each MCC (Fig. 1, black lines) is obtained from the well-known Chua's circuit [15], by replacing the nonlinear locally-active resistor (Chua's diode) with an ideal locallyactive flux-controlled memristor, whose constitutive relation is q = f (ϕ), where ϕ and q are the memristor flux and charge, respectively [16].
The state equations (SEs) in the (v, i)-domain of the MCCs array are (j = 1, . . . , N)  in the projection (v 1 , v 2 ) of the state space. Each plotted manifold admits multiple attractors: the top one admits a limit cycle (blue) and a chaotic attractor (yellow), the middle one two limit cycles (blue and yellow) and a chaotic attractor (red), the bottom one two limit cycles (blue and red).
where R jk denotes the resistor connecting the j − th and k − th MCCs through capacitors C 1,j and C 1,k . These resistors have the same value R jk = γρ, with γ = 1 if |j − k| = 1 and γ = 0 otherwise. We also assume boundary conditions of Dirichlet type, i.e., v 1, Note that (i) each invariant manifold is uniquely identified by [¯ 1 , . . . ,¯ N ] ∈ R N and (ii) each trajectory must lie on one and only one invariant manifold, as illustrated in Fig. 2. We remark that each invariant manifold is nonlinear, but in the considered projection of the state space it is a hyperplane. We then consider the MCC description in the (ϕ, q)-domain. By introducing the incremental flux and charge (j = 1, . . . , N) we obtain the following 3N coupled first-order SEs in the (ϕ, q)-domain Let us introduce the following normalized quantities: and let n(·) = Rf (·). By using the change of variables x j = ϕ j , y j = Li j , z j = Li j − ϕ j + RC 2 v 2,j , we obtain the dimensionless SEs of the considered array (j = 1, . . . , N): If we pick the initial conditions χ(0) for the state variables in the (v, i)-domain so that we are on the 0-manifold M(0), we have j0 = 0 and Eqs. (4) describe an array of memristor-less Chua's circuits with a nonlinear resistor having driving-point characteristic n(·). By contrast, if χ(0) is such that 0 = (χ(0)) = 0, the single MCC evolves on manifold M( 0 ) and Eqs. (4) describe the same array of memristor-less Chua's circuits, but with a non-zero constant forcing term α j0 , j = 1, . . . , N.

III. SYNCHRONIZATION
Our aim is to study the global synchronization of a 1D array of N diffusively-coupled identical MCCs in the (ϕ, q)-domain. We say that the network described by Eqs. (4) synchronizes The synchronization of the N systems shrinks the dimension of the subspace in which the synchronous solution lies, with respect to the dimension of the full state space. This subspace is called synchronization manifold and is defined as Since the coupling is diffusive, the necessary and sufficient condition for the invariance of M S [17] is 10 = . . . = N0 =¯ ∈ R. This means that the isolated MCCs are identical, as they have the same constant forcing term α j0 . To assess the stability of the synchronization manifold, i.e., to establish whether initial conditions x x x 1 (0) = . . . = x x x N (0) get attracted by the synchronous manifold or not, we rely on two different approaches, described in next sections, and on three related methods.

A. Analytic Conditions for Synchronization
We wish to apply the analytic results in [13] (in particular Corollary 3 and the analysis for a symmetric tridiagonal coupling) to ensure synchronization of the MCC array.
First, we have to find the minimum positive T 1 ∈ R such that the isolated MCC in the (ϕ, q)-domain with a negative feedback T = diag(T 1 , 0, 0) is uniformly asymptotically stable in the sense (Yoshizawa stability) defined in [13]: According to [13], T 1 simply needs to be not smaller than the absolute value of the largest negative slope of f (·). Therefore, a sufficient condition on T 1 is From formula (26) in [13], we obtain that a sufficient condition for the array of MCCs to globally synchronize is where σ is the normalized parameter introduced in Eq. Remark: The sufficient conditions on T 1 and σ are independent of j0 , namely, they ensure synchronization for any manifold and for any attractor on a given manifold. Actually, we expect that the conditions ensuring synchronization, which are related to the dynamics close to an attractor, depend upon both the manifold where a trajectory lies and also each attractor on a manifold, in the case of coexisting attractors; therefore, the above analytic conditions can be quite conservative.

B. Master Stability Function Approach
Note that plugging the synchronization condition x x x 1 Since the N systems have the same state and the same derivative, they will follow the same synchronized time evolution, described by the equatioṅ Eqs. (4) can also be rewritten more compactly aṡ where A = {a jk } is the Laplacian matrix, a jk = d jk −δ jk k d jk , and δ jk is the Kronecker's delta. The stability analysis of the synchronous solution through the Master Stability Function (MSF) approach [14], [18] is based on the study of small perturbations w w w j (τ ) = x x x j (τ )−x x x s (τ ) about the synchronous solution, each of which evolves based on the variational equation, for j = 1, . . . , N, where D is the Jacobian operator. Following the MSF approach, Eq. (11) can be decoupled into a number of independent blocks of the typė j = 1, . . . , N, where η η η j is the eigenmode associated with the eigenvalue λ j of A. Let us introduce the MSF M(ξ ), which associates the Maximum Lyapunov Exponent (MLE) of Eq. (12) to any possible value ξ of σ λ j . We remark that only the (N−1) perturbations η η η 2 (t), η η η 3 (t), . . . , η η η N (t) associated with the eigenvalues λ 2 , λ 3 , . . . , λ N determine the stability of the synchronization manifold. As a result, the synchronization manifold is stable if and only ifM(σ ) = max j M(σ λ j ) < 0 for j = 2, . . . , N.

A. Bifurcation Diagram for the Isolated MCC
In this section, we analyze the possible behaviors that the isolated MCC can exhibit, by setting α = 9.5 and β = 15. First, we define a regular grid of 40 points in the interval ∈ [0, 1]; then we compute the brute-force bifurcation diagram by simulating Eq. (4) starting from 20 different initial conditions. Depending on the considered manifold, for a given initial condition the isolated MCC can converge to one out of three different stable attractors (one periodic and two chaotic), as depicted in the Poincaré section in Fig. 3: for¯ < −0.75 MCC can converge to the blue (periodic) and red (chaotic) Fig. 3. Bifurcation diagram of an isolated MCC (intersections of asymptotic trajectories with a Poincaré section) vs.¯ . The circuit has up to three attractors (blue, red and yellow, with the same color code as in Fig. 2).

B. Global Synchronization in an Array of MCCs
We consider an array of N = 4 MCCs. To analyze the synchronization of this array, we compare three methods: Method A -This method is purely analytic. We obtain T 1 from Eq. (7) and the stability bound of the synchronous solution from Eq. (8).
Method B -This method is semi-analytic. We obtain T 1 by solving numerically Eq. (6) and the stability bound of the synchronous solution again from Eq. (8).
Method C -This method is purely numerical. We obtain the stability bound of the synchronous solution throughM(ξ ). Figure 4 shows the results of this analysis with respect to the invariant set (one per panel), the manifold (through¯ ), and the connection strength (through σ ). The functionM(ξ ) is color-coded as shown in the top bar: a shade of blue indicates the stability of the synchronous solution (i.e., negative MLE), whereas colors from green to yellow denote instability (i.e., positive MLE). The black line marks the stability edge obtained through method C. The dashed red vertical lines mark the sufficient condition obtained with method A. Notice that the condition is the same for any invariant set. The solid red lines mark the condition obtained with method B. Notice that, since the MSF grows monotonically with σ , the stability is determined only by the maximum eigenvalue ψ(N) of the connectivity matrix. For this reason, for N = 4 we get a similar diagram with a re-scaled x-axis, as shown in the bottom part of Fig. 4.

C. More Realistic Memristor Model
The proposed techniques can also be applied to more realistic memristor circuits. Here, we briefly discuss the relevant case where the ideal memristor is replaced by a real extended memristor [11] given by the parallel connection of an ideal flux-controlled memristor and a diode-like voltage-controlled nonlinear resistor i = h(v) [17] (see the gray part of Fig. 1), where the latter accounts for the rectification effect in the off-state due to the Schottky-like transport at one of the metal/oxide memristor interfaces [19]. The MCC with the new memristor model (called real-MCC) is shown in Fig. 1.
The first SE in the (v, i)-domain changes accordingly: where, to fix ideas, we let h(v) = I S (exp( v ηV T ) − 1), I S = 10 −12 A, V T = 26 · 10 −3 V and η = 1.7. The isolated oscillator is multistable, therefore its state can converge to an attractor that depends on the initial conditions. Notice that T 1 = a is still a stability boundary (method A) for the more realistic memristor model, even if largely conservative. This can be easily verified owing to the passivity of the added nonlinear resistor. We apply the proposed methods to an array of N = 4 real-MCCs; in particular, we analyze the stability of two attractors: a chaotic attractor (top panels of Fig. 5) and a 4-cycle, i.e., a limit cycle with four oscillations per period (bottom panels of Fig. 5). The right panels showM (blue curve) and the edge of stability computed with method A (dashed red line), method B (solid red), and method C (black) for the two considered attractors. We remark that this diagram holds for any N. It is apparent that also in this case method C provides the tightest condition, at the cost of a higher computational cost, as better discussed below.

V. CONCLUDING REMARKS
In this brief, we have provided a stability analysis of the synchronous solution in arrays of MCCs, by comparing three methods. Method A has a negligible computational cost (less than a second in the considered examples 1 ) but provides a conservative condition. It can be extended quite easily to other network topologies but does not allow investigating separately the synchronization of different coexisting attractors. Method B is more accurate, but computationally heavier (few seconds in the considered examples 1 ). Method C is the most accurate and is attractor-specific, but requires the heaviest computational effort and is network-dependent (few minutes in the considered examples 1 ). The obtained results could be extended from several standpoints. First of all, other kinds of synchronization (such as cluster synchronization [20], [21], [22], [23]) can be analyzed. Second, different network topologies can be investigated, either synthetic structures (such as all-to-all, circular array, and tridiagonal structure [13]) or structures related to real applications of memristors (regular lattice, nearest neighbor, and small world [24]). Last but not least, in the MCCs we considered both an ideal flux-controlled memristor, as in the original paper by Chua [16], and a real model accounting for interface physical effects. Other more realistic models for both the memristor [11] and the other network elements (including the coupling components) would require a higher computational effort for method C and further investigations to find analytic conditions for methods A and B. As a final remark, the proposed methods provide information about models of real-world circuits and systems and therefore qualitative information about their possible behaviors. To obtain more quantitative information, a further step is required, i.e., to obtain experimental bifurcation diagrams, thus switching from theory to practice [25], [26], [27], [28].