Emergence and suppression of chaos in small interdependent networks by localized periodic excitations

Abstract

Mastering the emergence and suppression of chaos of complex networks is currently a fundamental task for the nonlinear science community with potential relevant applications in diverse fields such as microgrid technologies, neural control engineering, and ecological networks. Here, the emergence and suppression of chaos in a complex network of driven damped pendula, which consists of two small starlike networks coupled by a single link, is investigated in the case where only the two hubs are subject to impulse-induced chaos-mastering excitations. Distinct chaos-mastering scenarios are found which depend on whether the connectivity strategy between the starlike networks is hub-to-hub, hub-to-leaf, or leaf-to-leaf. An explanation is given of the underlying physical mechanisms of these chaos-mastering scenarios and their main characteristics. The findings may be seen as a contribution to an intermediate step towards the long-term goal of mastering chaos in scale-free networks of damped-driven nonlinear systems.

1 Introduction

The ability to master complex networks of dissipative systems is of the uttermost relevance to many fundamental problems in science and technology [1,2,3,4,5,6], including dynamic pattern evolution [7], signal amplification [8], synchronization of individual dynamical behaviour occurring at a network’s vertices [9], and transport [10]. Particularly relevant because of its many potential applications is the domain of the order–chaos transition in the dynamics of a complex network, including neuronal networks [11] and microgrid technologies [12]. While complex networks are usually non-homogeneous structures, heterogeneous connectivity in the form of a scale-free topology [13] is a common property of many diverse real-world networks, which means that most nodes have few connections while just a small set of them are highly connected − the so-called hubs. Remarkably, starlike structures are the main motifs of scale-free networks, which has recently motivated the study of diverse synchronization phenomena in oscillator networks with starlike couplings [14, 15] as well as the suppression and emergence of chaos in starlike networks (SNs) of damped-driven nonlinear systems [16, 17]. Scale-free networks are ubiquitous examples of the general case where networks interact with other networks of similar or distinct natures, forming the so-called networks of networks (NONs) [18,19,20,21,22,23].

An essential property of complex dynamical networks is that perturbations to one node, for instance, in the form of periodic excitations, can modify the behaviour of other nodes, potentially entailing variation in the dynamics of the whole network. This means that local injection of energy can modify a network’s energy landscapes in parameter space and hence reshape the basins of attraction of the possible attractors in phase space. A similar thought may be exploited to master network dynamics [24] in distinct connected systems emerging in ecology [25], biology [26, 27], and physics [17, 28]. In particular, we shall study the manipulation of small-sized NONs in the sense of driving the NON from certain non-chaotic (chaotic) initial states to certain final stable chaotic (non-chaotic) states with the purpose of discovering not simply the driver nodes but also fixing their relative effectiveness as well as acquiring helpful information on the zones of parameter space in which appropriate signals in the way of parametric excitations (PEs) may be efficient. As an illuminating step towards the long-term aim of mastering chaos in scale-free networks, we investigate here the appearance and cancellation of chaos and conjoint synchronization phenomena by locally changing the impulse provided by PEs in the case of two connected small SNs as one of the simplest models of a NON. The impulse provided by a zero-mean periodic excitation F(t) having equidistant zeros (or excitation’s impulse for short,

$$\begin{aligned} I\equiv \left| \int _{0}^{T/2}F(t)dt\right| , \end{aligned}$$

with T being the period) is a quantity taking into account the combined effects of the excitation’s period, amplitude, and waveform. The value of the excitation’s impulse has been supported previously in quite diverse contexts, such as nonautonomous Hamiltonian systems [29], chaotic lasers [30], ratchets [31,32,33], discrete breathers in nonlinear oscillator networks [34], signal amplification in scale-free networks [35], and wave-packet localization in driven quantum systems [36], to quote just a few.

We aim here to explore a topology-induced chaos-order scenario in small two-hub NONs of damped-driven systems subjected to local chaos-taming PEs maintaining their common period and amplitude fixed. While impulse-induced emergence and suppression of chaos in isolated SNs has been investigated formerly in Ref.  [17], we shall here study such an impulse-induced chaos-order scenario in two-hub NONs for the instance of a single connector link (or interlink) connecting two SNs. In particular, the findings will be discussed through the analysis of two linked SNs of damped, parametrically excited pendula (DPEPs)—see Fig. 1a. The DPEP is a standard model that is sufficiently simple to allow analytical studies [17] while presenting universal characteristics, which are typical of numerous dissipative chaotic systems.

Fig. 1

a Schematic representations of three examples of small NONs, each of which consists of two SNs connected by a single interlink: hub-to-hub, hub-to-leaf, and leaf-to-leaf, respectively from left to right. b Parametric excitation function $f(t)\equiv A(m){\text {sn}}\left[ 4K(m)t/T;m\right] {\text {*}}{dn}\left[ 4K(m)t/T;m\right]$ [cf. Equation (2)] vs t/T for $m=0$ (thin line), $m=0.717\simeq m_{\max }^{impulse}$ (medium line), and $m=0.999$ (thick line). The inset shows the normalized impulse (see the text) vs the shape parameter. The quantities plotted are dimensionless

In this present communication, we explore the new chaos-order scenario emerging from the connection of two small SNs depicted by Eq. (1) in Sec. II by adopting sets of parameter values such that each isolated pendulum subjected to harmonic PEs in each SN exhibits non-chaotic behaviour. Next, we analyse the interplay between localized impulse-induced mastering of chaos and heterogeneous connectivity, which is due to the existence of an interlink and $N_{1}+N_{2}-2$ intralinks, when the shape parameter m is taken to be $m=0$ (i.e. harmonic PE) except for determined sets of pendula that are subjected to PEs of changeable waveform (i.e. $m\in \left] 0,1\right[$, cf. Equation (2)).

The rest of the communication is organized as follows: In Sec. II, we provide a description of the model equations and explore the interplay between local change of the PE’s impulse and heterogeneous connectivity in the simple instance of a NON composed of two SNs connected through a single interlink in accordance with three different connectivity schemes: hub-to-hub, hub-to-leaf, and leaf-to-leaf. We characterize the principal properties of a quite complex chaos-order scenario and establish how the effectiveness of the local reshaping of PEs at suppressing and inducing chaos depends upon the connectivity scheme of the two SNs, the degree of connectivity of the target nodes, and the values of the coupling constant associated with interlinks. Finally, Sec. III is devoted to a discussion of the major findings and of future work.

2 Localized impulse-induced chaos-taming in a network of networks

The model system of the two uncoupled SNs having $N_{1},N_{2}$ nodes, respectively, is:

$$\begin{aligned}&\overset{..}{\theta }_{H}+\left[ 1+\gamma f_{H}(t)\right] \sin \theta _{H} \nonumber \\&\quad =-\mu \overset{.}{\theta }_{H}+\lambda \sum _{i=1}^{N-1}\sin \left( \theta _{i}-\theta _{H}\right) ,\nonumber \\&\overset{..}{\theta }_{i}+\left[ 1+\gamma f_{i}(t)\right] \sin \theta _{i} \nonumber \\&\quad =-\mu \overset{.}{\theta }_{i}+\lambda \sin \left( \theta _{H}-\theta _{i}\right) , \end{aligned}$$

(1)

$i=1,...,N-1$, $N=\left\{ N_{1},N_{2}\right\}$, where $f_{H,i}(t)$ is a zero-mean unit-amplitude T-periodic PE:

$$\begin{aligned} f_{H,i}(t)&\equiv A(m){\text {sn}}\left( \varOmega _{H,i}t;m_{H,i} \right) {\text {*}}{dn}\left( \varOmega _{H,i}t;m_{H,i}\right) ,\nonumber \\ A(m)&\equiv \left[ \alpha _{1}+\frac{\alpha _{2}}{1+\exp \left[ \left( m-\alpha _{3}\right) /\alpha _{4}\right] }\right] ^{-1}, \end{aligned}$$

(2)

with $\alpha _{1}=0.43932,\alpha _{2}=0.69796,\alpha _{3}=0.37270,\alpha _{4}=0.26883$. Notice that all parameters and variables are dimensionless: $\varOmega _{H,i}=\varOmega _{H,i}\left( T,m_{H,i}\right) \equiv 4K(m_{H,i})/T$, $\gamma$ and T are the common excitation amplitude and period, respectively, $\lambda >0$ is the coupling constant, $\mu$ is the damping coefficient, ${\text {*}}{dn}\left( \cdot ;m\right) ,{\text {sn}}\left( \cdot ;m\right)$ are Jacobian elliptic functions of parameter m, A(m) is a normalization function which is introduced for the elliptic PE $f_{H,i}(t)$ to have the same period T and amplitude 1 for any waveform (i.e. $\forall m\in [ 0,1[$, see Fig. 1b), and K(m) is the complete elliptic integral of the first kind. It is worth mentioning that, contrary to other model systems such as the Kuramoto model, the coupling constant does not scale with N. Equations (1) depict the dynamics of a highly connected pendulum (or hub), $\theta _{H}$, and $N-1$ linked pendula (or leaves), $\theta _{i}$. The PE $f_{H,i}(t)$ was chosen for presenting three noteworthy features. First, when $m=0$ then $f_{H,i}(t)=\sin \left( 2\pi t/T\right)$, i.e. one retrieves the standard case of a sinusoidal PE [37,38,39,40,41], while, for the limiting value $m=1$, $f_{H,i}(t)$ vanishes. Second, its waveform (and hence its impulse) is varied by solely changing the shape parameter, m, between 0 and 1. And third, the PE’s impulse per unit of period, $I=A(m)/\left[ 2K(m)\right]$, exhibits a single maximum at a definite value of the shape parameter: $m=m_{\max }^{impulse}\simeq 0.717$.

We consider NONs consisting of two SNs of the type described by Eqs. (1)–(2), which are connected by a single interlink of the same type as that of the intralinks of the two SNs, but with different coupling constants: $\lambda _{HH},\lambda _{LL}$, and $\lambda _{HL}$ depending upon whether the strategy of connectivity between them is hub-to-hub, leaf-to-leaf, or hub-to-leaf, respectively, while $\lambda$ is the common coupling constant of the intralinks of the two connected SNs (see Fig. 1a). The dynamic equation of these NONs was integrated numerically using a fourth-order Runge–Kutta algorithm with time step $dt=T/100$. To visualize their overall spatiotemporal dynamics when the two SNs are coupled, we calculated the average velocity

$$\begin{aligned} \sigma \left( jT\right) \equiv \frac{1}{N_{1}+N_{2}}\sum _{i=1}^{N_{1}+N_{2} }\overset{.}{\theta }_{i}\left( jT\right) , \end{aligned}$$

(3)

where j is an integer multiple of the pulse period T, while their degree of synchronization is characterized by the correlation function [16]

$$\begin{aligned} C\equiv \frac{2}{\left( N_{1}+N_{2}\right) (N_{1}+N_{2}-1)}\sum _{\left( il\right) }\left\langle \cos \left( \theta _{i}-\theta _{l}\right) \right\rangle _{t},\nonumber \\ \end{aligned}$$

(4)

with the summation being over all pairs of pendula, and where $\left\langle \cdot \right\rangle _{t}$ indicates time averaging over a pre-defined (sufficiently long) observation window. Note that C is 1 for the perfectly synchronized state, while desynchronization increases as C decreases from 1. To check the reliability of C as a synchronization measure, we alternatively calculated the standard deviation

$$\begin{aligned} \varDelta \equiv \left\langle \sqrt{\frac{1}{N_{1}+N_{2}}\sum _{i=1}^{N_{1}+N_{2} }\left( \overset{.}{\theta }_{i}-\overline{\overset{.}{\theta }}\right) ^{2} }\right\rangle _{t} \end{aligned}$$

(5)

with $\overline{\overset{.}{\theta }}\equiv \left( N_{1}+N_{2}\right) ^{-1}\sum _{i=1}^{N_{1}+N_{2}}\overset{.}{\theta }_{i}$. We used bifurcation diagrams to illustrate our findings, constructing them by means of a Poincaré map at the parameters indicated in the figure captions. Starting at a certain value of the parameter of interest (for instance, $\lambda _{HH}=0$) and taking the transient time as $7.5\times 10^{4}$ PE periods after every increment (for instance, $\varDelta \lambda _{HH}=3.3\times 10^{-3}$), we sampled 30 PE periods by picking up the first $\sigma$ values of every pulse cycle, while to obtain the correlation function [cf. Equation (4)] and the standard deviation [cf. Equation (5)] we calculated C and $\varDelta$ averaged over 200 additional PE periods. The initial conditions were chosen randomly and independently for each node of the two coupled SNs, while the emergence and extent in parameter space of chaos (chaotic attractors) were characterized by using Lyapunov exponent (LE) calculations. We computed the LEs using a version of the algorithm introduced in Ref.  [42, 43] and typically integrated up to $7.5\times 10^{4}$ drive cycles for illustrative fixed parameters.

Before applying any chaos-taming excitation, we assume a set of fixed parameters $\left( \gamma ,\mu ,T\right)$ such that each isolated pendulum subjected to a sinusoidal PE $\left( m_{H}=m_{i}=0,i=1,...,N-1\right)$ displays regular behaviour. In particular, we consider the least favourable situation for the impulse-induced emergence of chaos: when the steady state of the isolated pendulum is the equilibrium $\left( \theta =0,\overset{.}{\theta }=0\right)$, i.e. the state of minimum energy. In other words, we study the impulse-induced destabilization of the synchronized equilibrium giving rise to generic chaotic states. Since we are mainly interested in studying the effectiveness of the localized impulse-induced injection of energy at giving rise to chaos in NONs, we shall focus in the present work on the least favourable situation for its emergence, i.e. when we are solely mastering the most highly connected pendulum in each SN (see Ref.  [17]). For the sake of clarity, we shall consider the three connectivity strategies separately.

2.1 Hub-to-hub interlink

The complete model of the NON reads:

$$\begin{aligned}&\overset{..}{\theta }_{H}^{(1)}+\left[ 1+\gamma f_{H}^{(1)}(t)\right] \sin \theta _{H}^{(1)} \nonumber \\&\quad =\lambda \sum _{i=1}^{N_{1}-1}\sin \left( \theta _{i}^{(1)}-\theta _{H}^{(1)}\right) \nonumber \\&\qquad +\lambda _{HH}\sin \left( \theta _{H}^{(2)}-\theta _{H}^{(1)}\right) -\mu \overset{.}{\theta }_{H}^{(1)},\nonumber \\&\overset{..}{\theta }_{i}^{(1)}+\left[ 1+\gamma f_{i}^{(1)}(t)\right] \sin \theta _{i}^{(1)} \nonumber \\&\quad =-\mu \overset{.}{\theta }_{i}^{(1)}+\lambda \sin \left( \theta _{H}^{(1)}-\theta _{i}^{(1)}\right) ,\nonumber \\&\overset{..}{\theta }_{H}^{(2)}+\left[ 1+\gamma f_{H}^{(2)}(t)\right] \sin \theta _{H}^{(2)} \nonumber \\&\quad =\lambda \sum _{j=1}^{N_{2}-1}\sin \left( \theta _{j}^{(2)}-\theta _{H}^{(2)}\right) \nonumber \\&\qquad +\lambda _{HH}\sin \left( \theta _{H}^{(1)}-\theta _{H}^{(2)}\right) -\mu \overset{.}{\theta }_{H}^{(2)},\nonumber \\&\overset{..}{\theta }_{j}^{(2)}+\left[ 1+\gamma f_{j}^{(2)}(t)\right] \sin \theta _{j}^{(2)}\nonumber \\&\quad =-\mu \overset{.}{\theta }_{j}^{(2)}+\lambda \sin \left( \theta _{H}^{(2)}-\theta _{j}^{(2)}\right) , \end{aligned}$$

(6)

$i=1,...,N_{1}-1$, $j=1,...,N_{2}-1$, where $f_{H}^{(1)}(t),f_{i} ^{(1)}(t),f_{H}^{(2)}(t),f_{i}^{(2)}(t)$ are the PEs given by Eq. (2) and $\lambda _{HH}>0$ is the coupling constant of the hub-to-hub interlink.

Let us consider first the simple symmetric case of two identical SNs, $N_{1}=N_{2}$. One finds that when the isolated SNs do not exhibit chaos for the same set of the remaining parameters, the same behaviour holds after coupling the two hubs for (almost) any value of the interlink coupling $\lambda _{HH}\in \left[ 0,1\right]$, as we verified for the cases $N_{1}=N_{2}=6$ and $N_{1}=N_{2}=5$ (see Figs. S6 and S7 in the Supplementary Information).

Fig. 2

Bifurcation diagrams of the average velocity $\sigma$ [cf. Equation (3)] and maximal Lyapunov exponent $\varLambda _{\max }$ (solid lines) as functions of the shape parameter $m=m_{1H}=m_{2H}$ for a NON with a hub-to-hub interlink [cf. Equation (6)], $\lambda _{HH}=0.01$ and three values of the numbers of pendula: a $N_{1}=4,N_{2}=6$, b $N_{1}=3,N_{2}=6$, and c $N_{1}=2,N_{2}=6$. Versions d, e, f show the respective corresponding bifurcation diagrams of the average velocity $\sigma$ and maximal Lyapunov exponent $\varLambda _{\max }$ as functions of the interlink coupling $\lambda _{HH}$ for $m=m_{1H} =m_{2H}=0.65$. Fixed parameters: $\lambda =0.01,\gamma =0.88,\mu =0.1,T=4.1$. The quantities plotted are dimensionless

The scenario is quite different when the hubs’ degrees $N_{1}-1,N_{2}-1$ are distinct, even in the situation where one of the SNs is in a regular or slightly chaotic state when it is isolated $\left( \lambda _{HH}=0\right)$, such as for the cases $N_{1}=4,N_{2}=6$, $N_{1}=3,N_{2}=6$, and $N_{1}=2,N_{2}=6$ [see Fig. 2a, b, and c, respectively]. We found that chaos emergence occurs over relatively small intervals of values of the shape parameter m, which are roughly centred around the value predicted from Melnikov method for an isolated pendulum $\left( m=m_{th}(T=4.1)\simeq 0.6562\right)$ [17], while the width of these chaotic ranges increases as the difference $\left| N_{1}-N_{2}\right|$ is increased [cf. Figure 2a, b and 2(c)]. One also sees that, as this difference is increased, the existence of windows of regularization interspersed with windows of chaos becomes more noticeable (compare Fig. 2a with Fig. 2c). Also, the chaotic behaviour, which is limited to the weak coupling regime for isolated SNs [17], occurs over wide ranges of the interlink coupling $\lambda _{HH}$ after linking the two SNs, such as for the cases $N_{1} =4,N_{2}=6$, $N_{1}=3,N_{2}=6$, and $N_{1}=2,N_{2}=6$ [see Fig. 2d, e and f, respectively]. Typically, the overall behaviour of the maximal LE, $\varLambda _{\max }$, as a function of $\lambda _{HH}$, presents a single maximum at a value $\lambda _{HH,\max }$ that is approximately proportional $\left( \lambda _{HH,\max }\sim 10\lambda \right)$ to the common value of the intralink coupling [$\lambda =0.01$ in Fig. 2d, e and f], almost independently of the degree of asymmetry of the NON $\left| N_{2}-N_{1}\right|$ (see Fig. S8 in the Supplementary Information). Thus, a particular coupling seems to favour the impulse-induced chaotic dynamics of the NON with hub-to-hub coupling. This enhancement effect is maximal when $m\approx m_{th}(T=4.1)\simeq 0.6562$, i.e. when the PE’s impulse is sufficiently high [see Figs. 2a, b and c], and occurs together with increasing desynchronization as the difference $\left| N_{1}-N_{2}\right|$ is increased [see Fig. 3a and c] and is again a consequence of breaking the symmetry of the NON.

Fig. 3

a Correlation function C [cf. Equation (4)] and c standard deviation $\varDelta$ [cf. Equation (5)] as functions of the shape parameter $m=m_{1H}=m_{2H}$ for a NON with a hub-to-hub interlink [cf. Equation (6)], $\lambda _{HH}=0.01$ and three values of the numbers of pendula: $N_{1}=4,N_{2}=6$ [thick red (grey) line]; $N_{1}=3,N_{2}=6$ [medium green (light grey) line]; $N_{1}=2,N_{2}=6$ [thin blue (black) line]. b Correlation function C and (d) standard deviation $\varDelta$ as functions of the interlink coupling $\lambda _{HH}$ for $m=m_{1H}=m_{2H}=0.65$ and three values of the numbers of pendula: $N_{1}=4,N_{2}=6$ [medium red (grey) line]; $N_{1}=3,N_{2}=6$ [thick green (light grey) line]; $N_{1}=2,N_{2}=6$ [thin blue (black) line]. Fixed parameters as in Fig. 2. The quantities plotted are dimensionless

Besides this chaos-inducing desynchronization around the value $\lambda _{HH}\approx \lambda _{HH,\max }\simeq \lambda$, Fig. 3b and d show that an even more drastic desynchronization occurs at greater values of the interlink coupling $\left( 0.7\lesssim \lambda _{HH}\lesssim 0.8\right)$, which are clearly beyond the weak coupling regime.

Fig. 4

Maximal Lyapunov exponent $\varLambda _{\max }$ (left) and correlation function C [cf. Equation (4)] (right) in the $m-\lambda _{HH}$ parameter plane for $N_{1}=4,N_{2}=6$ (top panel), $N_{1}=3,N_{2}=6$ (middle panel), and $N_{1}=2,N_{2}=6$ (bottom panel), respectively. Fixed parameters as in Fig. 2. The quantities plotted are dimensionless

Figure 4 left and right shows that the range of the shape parameter m roughly centred around $m=m_{th}(T=4.1)\simeq 0.6562$, which is where chaos and chaos-induced desynchronization appear, shrinks rapidly as the interlink coupling $\lambda _{HH}$ grows from zero, its width reaching a constant value over a relatively wide range. It is observed that the drastic desynchronization occurring in the range $0.7\lesssim \lambda _{HH}\lesssim 0.8$ is associated with a decrease in the maximal LE (cf. Figure 4 left). Worthy of note is the existence of additional desynchronization ranges that appear when the impulse transmitted by the PEs is sufficiently small (i.e. when m is close to 1, cf. Figure 4 right), irrespective of the values of $\lambda _{HH}$, while they are not associated with chaotic behaviour (cf. Figure 4 left). Figure 5 shows the angular velocity of all coupled pendula versus time after $4.5\times 10^{4}$ driving cycles for several values of the interlink coupling, shape parameter, and number of nodes of the two coupled SNs. In general, we found that desynchronization states appear due to cluster synchronization of different sets of pendula, including the two hubs [cf. Figure 5c, f and i], all or a number of the leaves of an SN [cf. Figure 5b and c], and all the leaves of the two SNs [cf. Figure 5f and i]. In the weak coupling regime, such as for the case $\lambda _{HH}=0.01$[cf. Figure 5a, b, d, e, g, and h], the NON typically presents complete (global) synchronization when the shape parameter is sufficiently close to 0 (i.e. when the PE’s waveform is sufficiently close to the harmonic waveform, see Fig. 5j), and the degree of asymmetry of the NON is sufficiently small, as in the case shown in Fig. 5e for $N_{1}=4,N_{2}=6$. This global synchronization breaks after increasing the degree of asymmetry $\left| N_{2} -N_{1}\right|$, such as for the case $N_{1}=3,N_{2}=6$ [cf. Figure 5d] where the two sets of leaves present synchronized (but distinct) intermittent states, while the two hubs present desynchronized chaotic states. A further increase in the degree of asymmetry $\left| N_{2}-N_{1}\right|$ results in the almost complete chaotification of the NON, such as for the case $N_{1}=2,N_{2}=6$ [cf. Figure 5a] where the hub of the smallest SN presents a regular (period-8) state and the leaves of the two SNs present phase-locked chaos. Beyond the weak coupling regime, one may find, for certain values of the interlink coupling, synchronization of the chaotic states of the two hubs, while some leaves change from intermittent behaviour to chaos as $\left| N_{2}-N_{1}\right|$ is increased, such as for the case $\lambda _{HH}=0.791$ [cf. Figure 5c, f, and i]. Some illustrative examples of the dependence of the maximal LE upon the parameters $T,\gamma ,\mu$ for several NONs are given in Figs. S1–S5 (Supplementary Information). In particular, one sees that by increasing the degree of asymmetry of the NON $\left( \left| N_{2}-N_{1}\right| \right)$, higher values of the amplitude are needed to achieve chaotic behaviour (see Fig. S4), on the one hand, and chaotic behaviour is maintained for higher values of the dissipation coefficient (see Fig. S5), on the other hand.

Fig. 5

Angular velocities $\dot{\theta }_{i}$ of the coupled pendula vs t/T [cf. Equation (6)] after $4.5\times 10^{3}$ driving cycles for a, b, c $N_{1}=2,N_{2}=6$, d, e, f $N_{1}=3,N_{2}=6$, and (g), (h), (i) $N_{1}=4,N_{2}=6$. a, d, g: $m=0.49,\lambda _{HH}=0.01$. b, e, h: $m=0.6,\lambda _{HH}=0.01$. c, f, i: $m=0.65,\lambda _{HH}=0.791$. $H_{i}$ and $L_{ij}$ denote the hub and the leaf j of SN i, respectively. j Parametric excitation function $f(t)\equiv A(m){\text {sn}}\left[ 4K(m)t/T;m\right] {\text {*}}{dn}\left[ 4K(m)t/T;m\right]$ [cf. Equation (2)] vs t/T for $m=0$ (dotted line), $m=0.49$ (dashed line), and $m=0.6$ (solid line). Fixed parameters as in Fig. 2. The quantities plotted are dimensionless

Physically, the above findings indicate that a relatively small chaotic SN is not only behaving as an effective energy source for the other SN when the PE’s impulse is sufficiently large to induce chaos, but its hub and the hub of the other SN together form an energy dimer-source, which is able to impose in general (but not in all cases) its chaotic behaviour on the remaining pendula of the NON over almost the entire range of values of the interlink coupling. This provides the microscopic mechanism leading to chaos in the NON described by Eq. (6). Note that this symmetry-breaking interlink-induced emergence of chaos is absent for certain particular values of $\lambda _{HH}$ when the NON’s symmetry is sufficiently broken [compare the cases $N_{1}=4,N_{2}=6$ and $N_{1}=2,N_{2}=6$; cf. Figures 2d and f, respectively].

2.2 Hub-to-leaf interlink

In this case, the complete model of the NON reads:

$$\begin{aligned}&\overset{..}{\theta }_{H}^{(1)}+\left[ 1+\gamma f_{H}^{(1)}(t)\right] \sin \theta _{H}^{(1)} \nonumber \\&\quad =\lambda \sum _{i=1}^{N_{1}-1}\sin \left( \theta _{i}^{(1)}-\theta _{H}^{(1)}\right) \nonumber \\&\qquad +\lambda _{HL}\sin \left( \theta _{k}^{(2)}-\theta _{H}^{(1)}\right) -\mu \overset{.}{\theta }_{H}^{(1)},\nonumber \\&\overset{..}{\theta }_{i}^{(1)}+\left[ 1+\gamma f_{i}^{(1)}(t)\right] \sin \theta _{i}^{(1)} \nonumber \\&\quad =-\mu \overset{.}{\theta }_{i}^{(1)}+\lambda \sin \left( \theta _{H}^{(1)}-\theta _{i}^{(1)}\right) ,\nonumber \\&\overset{..}{\theta }_{H}^{(2)}+\left[ 1+\gamma f_{H}^{(2)}(t)\right] \sin \theta _{H}^{(2)} \nonumber \\&\quad =\lambda \sum _{j=1}^{N_{2}-1}\sin \left( \theta _{j}^{(2)}-\theta _{H}^{(2)}\right) -\mu \overset{.}{\theta }_{H}^{(2)},\nonumber \\&\overset{..}{\theta }_{j}^{(2)}+\left[ 1+\gamma f_{j}^{(2)}(t)\right] \sin \theta _{j}^{(2)} \nonumber \\&\quad =-\mu \overset{.}{\theta }_{j}^{(2)}+\lambda \sin \left( \theta _{H}^{(2)}-\theta _{j}^{(2)}\right) \nonumber \\&\qquad +\delta _{jk}\lambda _{HL}\sin \left( \theta _{H}^{(1)}-\theta _{k} ^{(2)}\right) \end{aligned}$$

(7)

Fig. 6

Bifurcation diagrams of the average velocity $\sigma$ [cf. Equation (3)] and maximal Lyapunov exponent $\varLambda _{\max }$ (solid lines) as functions of the shape parameter $m=m_{1H}=m_{2H}$ for a NON with a hub-to-leaf interlink [cf. Equation (6)], $\lambda _{HL}=0.01$ and three values of the numbers of pendula: a $N_{1}=4,N_{2}=6$, b $N_{1}=3,N_{2}=6$, and c $N_{1}=2,N_{2}=6$. Versions d–f show the respective corresponding bifurcation diagrams of the average velocity $\sigma$ and maximal Lyapunov exponent $\varLambda _{\max }$ as functions of the interlink coupling $\lambda _{HL}$ for $m=m_{1H} =m_{2H}=0.65$. Fixed parameters as in Fig. 2. The quantities plotted are dimensionless

$i=1,...,N_{1}-1$, $j=1,...,N_{2}-1$, where k denotes the single peripheral pendulum of SN 2, which is connected to the central pendulum of SN 1 $\left( 1\leqslant k\leqslant N_{2}-1\right)$, $\delta _{nm}$ is the Kronecker delta, and $\lambda _{HL}>0$ is the coupling constant of the hub-to-leaf interlink.

Fig. 7

a Correlation function C [cf. Equation (4)] and c standard deviation $\varDelta$ [cf. Equation (5)] as functions of the shape parameter $m=m_{1H}=m_{2H}$ for a NON with a leaf-to-leaf interlink [cf. Equation (8)], $\lambda _{LL}=0.01$ and three values of the numbers of pendula: $N_{1}=4,N_{2}=6$ (thick line); $N_{1}=3,N_{2}=6$ (medium line); $N_{1}=2,N_{2}=6$ (thin line). b Correlation function C and d standard deviation $\varDelta$ as functions of the interlink coupling $\lambda _{LL}$ for $m=m_{1H}=m_{2H}=0.65$ and three values of the numbers of pendula: $N_{1}=4,N_{2}=6$ (thick line); $N_{1}=3,N_{2}=6$ (medium line); $N_{1}=2,N_{2}=6$ (thin line). Fixed parameters as in Fig. 2. The quantities plotted are dimensionless

Fig. 8

Maximal Lyapunov exponent $\varLambda _{\max }$ (left) and correlation function C [cf. Equation (4)] (right) in the $m-\lambda _{HL}$ parameter plane for $N_{1}=4,N_{2}=6$ (top panel), $N_{1}=3,N_{2}=6$ (middle panel), and $N_{1}=2,N_{2}=6$ (bottom panel), respectively. Fixed parameters as in Fig. 2. The quantities plotted are dimensionless

Three main features differentiate the present case from the above scenario for a hub-to-hub interlink.

First, in the simple symmetric case of two identical SNs, $N_{1}=N_{2}$, the hub-to-leaf interlink breaks the NON’s symmetry since the peripheral pendulum of SN 2, which is connected to the central pendulum of SN 1, may also be effectively considered to be a peripheral pendulum of the latter after connecting the SNs, thereby effectively increasing its degree by one unit, $N_{1}-1\rightarrow N_{1}$, while $N_{2}$ remains unchanged [see Fig. 1a and 6]. Indeed, one finds that when the (effective) isolated SNs do not exhibit chaos for the same set of the remaining parameters [17], the same behaviour holds after coupling them for almost any value of the interlink coupling $\lambda _{HL}\in \left[ 0,1\right]$, as we verified for the case $N_{1}=5,N_{2}=6$ (see Fig. S9 in the Supplementary Information).

Second, the chaos-transmitting effectiveness of the hub-to-leaf interlink is greater than that of the hub-to-hub interlink in the weak homogeneous coupling regime $\left( 0<\lambda _{HL}\approx \lambda \ll 1\right)$ due to the absence of an energy dimer-source. Instead, there exist two energy monomer sources connected by the H-L-H pathway, thus forming a more powerful energy chain-source, as for the case $N_{1}=3,N_{2}=6$ [see Fig. 6b] which is in contrast to the (equivalent) hub-to-hub interlink scenario case $N_{1}=4,N_{2}=6$ [see Fig. 2a]. Indeed, there is a wider range of the shape parameter for chaos in the former case than in the latter case. Remarkably, the strength of this enhancement effect is maximal when $m\approx m_{th}(T=4.1)\simeq 0.6562$, as indicated by the dependence of the maximal LE $\varLambda _{\max }$ on the shape parameter m [cf. Figures 6b and 2a], and occurs together with increasing desynchronization as the difference $\left| N_{1}-N_{2}\right|$ is increased [see Fig. 7a and c].

Third, the chaotic behaviour, which remains limited to the weak coupling regime of the interlink coupling $\lambda _{HL}$ when the degree of asymmetry of the NON $\left| N_{2}-N_{1}\right|$ is sufficiently small, undergoes a remarkable extension to greater values of $\lambda _{HL}$, and occurs together with increasing desynchronization as $\left| N_{2} -N_{1}\right|$ is increased [cf. Figure 6d, e, f, 7b, and 7d].

Figure 8 left and right shows that the range of the shape parameter m roughly centred around $m=m_{th}(T=4.1)\simeq 0.6562$, which is where chaos and chaos-induced desynchronization appear, shrinks rapidly as the interlink coupling $\lambda _{HL}$ grows from zero, becoming zero in width from a certain value of the interlink coupling. This is in contrast to the behaviour observed for the case of a hub-to-hub interlink (cf. Sec. II A). Note also the existence of additional desynchronization ranges that appear over the interval $m\in \left[ 0.33,0.43\right]$ for certain values of the interlink coupling belonging to the interval $\lambda _{HL}\in \left[ 0.3,0.4\right]$ (cf. Figure 8 right), while they are associated with non-chaotic dynamics (cf. Figure 8 left). Similarly to the case of a hub-to-hub interlink, desynchronization states typically appear due to cluster synchronization of different sets of pendula, including mainly all or a number of the leaves of an SN [cf. Figure 9a, b, d, e, and f], while the NON presents complete synchronization when the shape parameter is sufficiently far from $m\approx m_{th}(T=4.1)\simeq 0.6562$ and the degree of asymmetry of the NON is sufficiently small, as in the case shown in Fig. 9c for $N_{1}=4,N_{2}=6$.

Fig. 9

Angular velocities $\dot{\theta }_{i}$ of the coupled pendula vs t/T [cf. Equation (7)] after $4.5\times 10^{3}$ driving cycles for a, d $N_{1}=2,N_{2}=6$, b, e $N_{1}=3,N_{2}=6$, and c, f $N_{1}=4,N_{2}=6$. a, b, c: $m=0.7,\lambda _{HL}=0.02$. d, e, f: $m=0.63,\lambda _{HL}=0.01$. $H_{i}$ and $L_{ij}$ denote the hub and the leaf j of SN i, respectively. Fixed parameters as in Fig. 2. The hub $H_{1}$ is connected to the leaf $L_{21}$. The quantities plotted are dimensionless

2.3 Leaf-to-leaf interlink

In this remaining case, the complete model of the NON reads

$$\begin{aligned}&\overset{..}{\theta }_{H}^{(1)}+\left[ 1+\gamma f_{H}^{(1)}(t)\right] \sin \theta _{H}^{(1)} \nonumber \\&\quad =\lambda \sum _{i=1}^{N_{1}-1}\sin \left( \theta _{i}^{(1)}-\theta _{H}^{(1)}\right) -\mu \overset{.}{\theta }_{H}^{(1)},\nonumber \\&\overset{..}{\theta }_{i}^{(1)}+\left[ 1+\gamma f_{i}^{(1)}(t)\right] \sin \theta _{i}^{(1)} \nonumber \\&\quad =-\mu \overset{.}{\theta }_{i}^{(1)}+\lambda \sin \left( \theta _{H}^{(1)}-\theta _{i}^{(1)}\right) \nonumber \\&\qquad +\delta _{il}\lambda _{LL}\sin \left( \theta _{k}^{(2)}-\theta _{l} ^{(1)}\right) ,\nonumber \\&\overset{..}{\theta }_{H}^{(2)}+\left[ 1+\gamma f_{H}^{(2)}(t)\right] \sin \theta _{H}^{(2)} \nonumber \\&\quad =+\lambda \sum _{j=1}^{N_{2}-1}\sin \left( \theta _{j}^{(2)}-\theta _{H}^{(2)}\right) -\mu \overset{.}{\theta }_{H}^{(2)},\nonumber \\&\overset{..}{\theta }_{j}^{(2)}+\left[ 1+\gamma f_{j}^{(2)}(t)\right] \sin \theta _{j}^{(2)} \nonumber \\&\quad =-\mu \overset{.}{\theta }_{j}^{(2)}+\lambda \sin \left( \theta _{H}^{(2)}-\theta _{j}^{(2)}\right) \nonumber \\&\qquad +\delta _{jk}\lambda _{LL}\sin \left( \theta _{l}^{(1)}-\theta _{k} ^{(2)}\right) , \end{aligned}$$

(8)

$i=1,...,N_{1}-1$, $j=1,...,N_{2}-1$, where l, k denote the peripheral pendula connecting the two SNs ($1\leqslant k\leqslant N_{2}-1$, $1\leqslant l\leqslant N_{1}-1$), $\delta _{nm}$ is the Kronecker delta, and $\lambda _{LL}>0$ is the coupling constant of the leaf-to-leaf interlink.

Fig. 10

Bifurcation diagrams of the average velocity $\sigma$ [cf. Equation (3)] and maximal Lyapunov exponent $\varLambda _{\max }$ (solid lines) as functions of the shape parameter $m=m_{1H}=m_{2H}$ for a NON with a leaf-to-leaf interlink [cf. Equation (7)], $\lambda _{LL}=0.01$ and three values of the numbers of pendula: (a) $N_{1}=4,N_{2}=6$, (b) $N_{1}=3,N_{2}=6$, and (c) $N_{1}=2,N_{2}=6$. Versions (d), (e), (f) show the respective corresponding bifurcation diagrams of the average velocity $\sigma$ and maximal Lyapunov exponent $\varLambda _{\max }$ (solid lines) as functions of the interlink coupling $\lambda _{LL}$ for $m=m_{1H}=m_{2H}=0.65$. Fixed parameters as in Fig. 2. The quantities plotted are dimensionless

Let us consider first the simple symmetric case of two identical SNs, $N_{1}=N_{2}$. Similarly to the same case for a hub-to-hub interlink, one finds that when the isolated SNs do not exhibit chaos for the same set of the remaining parameters, the same behaviour holds after coupling the two leaves for almost any value of the interlink coupling $\lambda _{LL}\in \left[ 0,1\right]$, such as for the case $N_{1}=6,N_{2}=6$ (see Figs. S10 and S11 in the Supplementary Information). The scenario is quite different when the hubs’ degrees $N_{1},N_{2}$ are distinct. Thus, the chaos-transmitting effectiveness of the leaf-to-leaf interlink is greater than that of either the hub-to-hub or the hub-to-leaf interlinks due to the existence of two energy monomer sources connected by the H-L-L-H pathway, forming thereby a longer energy chain-source than in the other two types of interlinks, such as for the cases $N_{1}=4,N_{2}=6$, $N_{1}=3,N_{2}=6$, and $N_{1} =2,N_{2}=6$ [compare Fig. 2a, b, and c with Fig. 10a, b, and c, respectively]. Note, however, that the number and width of regularization windows interspersed within the chaotic interval roughly centred around $m=m_{th}(T=4.1)\simeq 0.6562$ are greater in the leaf-to-leaf interlink scenario than in the other two scenarios [cf. Figure 2c and c]. Similarly to the case of the hub-to-hub interlink, the chaotic behaviour occurs over wide ranges of the interlink coupling $\lambda _{LL}$ after linking the two SNs, such as for the cases $N_{1}=4,N_{2}=6$, $N_{1}=3,N_{2}=6$, and $N_{1}=2,N_{2}=6$ [see Fig. 10d, e, and f, respectively], and occurs again with increasing desynchronization as the difference $\left| N_{1}-N_{2}\right|$ is increased [see Fig. 11a and c]. Such a dependence of the chaos on $\lambda _{LL}$ is noticeably more uniform and of lesser strength than in the case of the hub-to-hub interlink (compare Fig. 4 left and 12 left), in the sense that the overall behaviour of the maximal LE $\varLambda _{\max }$ on average remains constant as a function of $\lambda _{LL}$, with the corresponding chaotic orbits being ever more uniform as the degree of asymmetry of the NON $\left| N_{2}-N_{1}\right|$ is increased [compare Fig. 10d, e, and f]. This uniformity is also reflected in the dependence of the correlation function C and the standard deviation $\varDelta$ on the interlink coupling $\lambda _{LL}$ [see Fig. 11b, d, and 12 right].

Fig. 11

a Correlation function C [cf. Equation (4)] and c standard deviation $\varDelta$ [cf. Equation (5)] as functions of the shape parameter $m=m_{1H}=m_{2H}$ for a NON with a leaf-to-leaf interlink [cf. Equation (8)], $\lambda _{LL}=0.01$ and three values of the numbers of pendula: $N_{1}=4,N_{2}=6$ (thick line); $N_{1}=3,N_{2}=6$ (medium line); $N_{1}=2,N_{2}=6$ (thin line). (b) Correlation function C and (d) standard deviation $\varDelta$ as functions of the interlink coupling $\lambda _{LL}$ for $m=m_{1H}=m_{2H}=0.65$ and three values of the numbers of pendula: $N_{1}=4,N_{2}=6$ (thick line); $N_{1}=3,N_{2}=6$ (medium line); $N_{1}=2,N_{2}=6$ (thin line). Fixed parameters as in Fig. 2. The quantities plotted are dimensionless

Also, we found the existence of additional desynchronization ranges that appear over ever larger intervals of the shape parameter as the degree of asymmetry of the NON $\left| N_{2}-N_{1}\right|$ is decreased and solely for certain values of the interlink coupling over the range $\lambda _{LL}\in \left[ 0.35,0.5\right] \cup \left[ 0.75,0.8\right]$ (cf. Figure 12 right), while they are not associated with chaos (cf. Figure 12 left). Similarly to the cases of the other two types of interlinks (cf. Secs. II A and II B), desynchronization states typically appear again due to cluster synchronization of different sets of pendula, including all or a number of the leaves of an SN [cf. Figure 13a, b, and e] and a number of the leaves of the two SNs [cf. Figure 13d and f], while the NON presents complete synchronization when the shape parameter is sufficiently far from $m\approx m_{th}(T=4.1)\simeq 0.6562$ and the degree of asymmetry of the NON is sufficiently small, as in the case shown in Fig. 13c for $N_{1}=4,N_{2}=6$.

Fig. 12

Maximal Lyapunov exponent $\varLambda _{\max }$ (left) and correlation function C [cf. Equation (4)] (right) in the $m-\lambda _{LL}$ parameter plane for $N_{1}=4,N_{2}=6$ (top panel), $N_{1}=3,N_{2}=6$ (middle panel), and $N_{1}=2,N_{2}=6$ (bottom panel), respectively. Fixed parameters as in Fig. 2. The quantities plotted are dimensionless

Fig. 13

Angular velocities $\dot{\theta }_{i}$ of the coupled pendula vs t/T [cf. Equation (8)] after $4.5\times 10^{3}$ driving cycles for (a), (b) $N_{1}=2,N_{2}=6$, (c), (d) $N_{1}=3,N_{2}=6$, and (e), (f) $N_{1}=4,N_{2}=6$. a, c, e: $m=0.49,\lambda _{LL}=0.02$. b, d, f: $m=0.65,\lambda _{LL}=0.09$. $H_{i}$ and $L_{ij}$ denote the hub and the leaf j of SN i, respectively. Fixed parameters as in Fig. 2. The leaf $L_{21}$ is connected to the leaf $L_{11}$. The quantities plotted are dimensionless

Physically, the aforementioned distinct effectiveness of the different interlink scenarios seems to stem from the distinct geodesic distance between the two target nodes (hubs): whereas the existence of an energy dimer-source corresponds to a minimal geodesic distance between the hubs (case of a hub-to-hub interlink), the existence of two energy monomer sources connected by the H-L-L-H pathway corresponds to a maximal distance between them (case of a leaf-to-leaf interlink).

3 Conclusions

In conclusion, we have shown how effective locally mastering the impulse transmitted by parametric periodic excitations is at inducing and taming chaos in two single-interlink-connected starlike networks of damped-driven pendula, leading to asynchronous chaotic states and equilibria, respectively. Three different chaos-mastering scenarios have been characterized depending upon whether the connectivity strategy between the starlike networks is hub-to-hub, hub-to-leaf, or leaf-to-leaf in the least favourable case for the emergence of chaos in which only the two hubs are subjected to impulse-induced chaos-mastering excitations, while each isolated pendulum subjected to a sinusoidal parametric excitation presents regular behaviour. When the parametric excitation’s impulse is large enough to induce chaos, the two pendula subjected to chaos-mastering excitations behave like energy sources for the remaining pendula, and the interplay between these two energy sources represents the microscopic physical mechanism leading to the whole network’s chaotification. We found that the distinct chaos-transmitting effectiveness of the different interlinks stems from the distinct geodesic distance between the two hubs (target nodes). Specifically, an energy dimer-source corresponds to a minimal geodesic distance between the target nodes (case of a hub-to-hub interlink), two energy monomer sources connected by the H-L-L-H pathway corresponds to a maximal distance between them (case of a leaf-to-leaf interlink), while the case of a hub-to-leaf interlink is associated with an intermediate distance between the connected hubs by the H-L-H pathway. Remarkably, we found that the emergence of chaos systematically occurs over ranges of the shape parameter roughly centred on the value predicted from the Melnikov method for an isolated pendulum [17], while these chaotic ranges typically contain the value of the shape parameter at which the impulse transmitted by the parametric excitation is maximum over wide ranges of the interlink coupling, irrespective of the connectivity strategy between the two starlike networks, thus providing additional evidence for the universal character of the impulse effect. The width of these chaotic ranges systematically increases as the degree of asymmetry of the whole network is increased, the chaos-inducing desynchronization of the whole network being typically maximum in the weak coupling regime and appearing due to cluster synchronization of different sets of pendula.

Finally, we hope that our findings can serve as a significant step towards achieving reliable mastering of chaos in larger complex interconnected networks. Since a highly connected node in a scale-free network can be thought of as a hub of a locally starlike part of the network, and the aforementioned microscopic physical mechanisms are expected to remain valid, we are presently exploring the effectiveness of the impulse-induced mastering of chaos in such complex networks. Another problem that deserves to be investigated is the robustness of the present impulse-induced chaos-mastering scenario against the presence of time-varying coupling in temporal networks [44].