Amplitude modulation control of spatiotemporal chaos in starlike networks of damped-driven pendula

Physica D 457 (2024) 133950 Available online 11 October 2023 0167-2789/© 2024 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Physica D journal homepage: www.elsevier.com/locate/physd Amplitude modulation control of spatiotemporal chaos in starlike networks of damped-driven pendula R. Chacón a,b, A. Martínez García-Hozc, F. Palmero d,∗ aDepartamento de Física Aplicada, Escuela de Ingenierías Industriales, Universidad de Extremadura, Apartado Postal 382, E-06006 Badajoz, Spain bInstituto de Computación Científica Avanzada (ICCAEx), Universidad de Extremadura, E-06006 Badajoz, Spain cDepartamento de Física Aplicada, E. T. S. de Ingenieros Agrónomos, Universidad de Castilla-La Mancha, E-13071 Ciudad Real, Spain dGrupo de Física No Lineal, Departamento de Física Aplicada I, Escuela Técnica Superior de Ingeniería Informática, Universidad de Sevilla, Avda Reina Mercedes s/n, E-41012 Sevilla, Spain ARTICLE INFO Communicated by T. Insperger Keywords: Chaos Control of chaos Parametrically excited damped pendulum Starlike networks of sinusoidally coupled chaotic pendula ABSTRACT Applying amplitude modulations to a parametrically excited damped pendulum is shown to be a reliable method to control (suppress or enhance) its chaotic behaviour. Analytical (Melnikov analysis) and numerical (Lyapunov exponents and bifurcation diagrams) results show an effective control scenario for a wide range of resonances between the two excitations implicated. Different routes of regularization as the chaos-controlling parameters vary are identified, including period-doubling and crises. The method’s effectiveness at suppressing spatiotemporal chaos of starlike networks of sinusoidally coupled chaotic pendula is demonstrated where effective regularization is obtained under localized control on an increasing number of pendula. 1. Introduction Since the end of the last century, the suppression and enhancement of the chaos emerging from general damped-driven systems have been a fundamental field of intensive investigation in a large number of contexts in science and engineering. Proof of this is the plentiful literature on the subject [1–4], including experimental research on chaos suppression [5–11]. The effectiveness of diverse suppressory techniques has been explored, including nonfeedback methods using small periodic excitations [12–16]. These latter methods have proven to be useful in applications where feedback control techniques are either inadequate or impractical, while they are easily applicable experimentally [10,16]. Furthermore, since three decades ago or so there has been the growth of a theoretical proposal to tame chaotic behaviour in dampeddriven systems which are capable of being studied by Melnikov’s method (MM), and which entails applying periodic chaos-controlling (CC) excitations [3,17]. This MM-based proposal has been demonstrated to be efficient in suppressing chaos in networks of damped, periodically driven, nonlinear oscillators [18,19]. Hitherto, mainly three kinds of instances have been explored regarding the two additive periodic excitations implicated (one chaos-inducing, CI, and the other CC): (i) both perturbations are forcings (external excitations), (ii) one excitation is parametric while the other is a forcing, and (iii) both perturbations are parametric excitations. To the best of our knowledge, however, the physically pertinent instance in which the CC ∗Corresponding author. E-mail address: [email protected] (F. Palmero). excitation acts multiplying the CI excitation has yet to be investigated in detail, despite its great theoretical relevance. Indeed, the application of small subharmonic amplitude modulations has been shown to be an effective procedure for preserving chaos in the context of laser systems [15,16], while the effectiveness of resonant amplitude modulations at suppressing chaos has been demonstrated for different damped nonlinear oscillators subjected to a sinusoidal forcing in the case where the two periodic excitations implicated are in phase [20–22]. Also, the effectiveness of a multiplicative method at controlling escape from a potential well by periodically modulating the amplitude of the escape-inducing excitations has been shown [23]. While it has been demonstrated that the main resonance between the two excitations implicated is the most effective resonance for optimally controlling homoclinic chaos in general [3], we shall show here that multiplicative CC excitations can tame chaotic behaviour beyond the main resonance, thus providing an additional technological motivation for the consideration of CC multiplicative excitations. Indeed, they could be useful for example in the problem of extracting energy from sea waves by means of a parametric pendulum [24]. In the present communication, we investigate this multiplicative CC scenario numerically and analytically by focusing on the paradigmatic model of a dissipative pendulum subjected to a parametric harmonic excitation [25–29]. Specifically, we demonstrate the effectiveness of the multiplicative control method in suppressing spatiotemporal chaos https://doi.org/10.1016/j.physd.2023.133950 Received 23 May 2023; Received in revised form 29 July 2023; Accepted 6 October 2023 Physica D: Nonlinear Phenomena 457 (2024) 133950 2 R. Chacón et al. of starlike networks of identical sinusoidally coupled chaotic pendula by localized application of CC excitations. It is worth mentioning that starlike networks are network motifs that occur repeatedly in real complex networks such as scale-free networks [30]. The remainder of the communication is organized as follows. Section 2discusses the MM- based analytical predictions for the dissipative parametrically excited pendulum subjected to small resonant amplitude modulations, and then contrasts numerical outcomes (bifurcation diagrams and Lyapunov exponents) with the analytical predictions, as well as describing the regularization routes in the control parameter plane for the case of resonances between the CI and CC excitations. Section 3explores the efficiency of the multiplicative control method under localized application of the CC excitations in starlike networks of sinusoidally coupled chaotic pendula. We consider the application of the CC excitations onto an increasing number of pendula, and study the robustness of the control scenario and related synchronization phenomena under changes in the number and degree of the controlled nodes as well as in the coupling strength. Finally, Section 4is devoted to a discussion of the major findings and of future work. 2. Multiplicative control scenario We consider the equation of a dissipative pendulum whose pivot is subjected to a vertical oscillation having a small resonant amplitude modulation: .. 𝜃+ sin 𝜃= −𝛿 . 𝜃−𝛾[1 + 𝜀cos (𝛺𝑡 +𝜙)]cos (𝜔𝑡)sin 𝜃, (1) where all variables and parameters are dimensionless. Here, the term describing the vertical sinusoidal oscillations of the pivot, −𝛾cos (𝜔𝜏) sin 𝜃, is to be regarded as the CI excitation, while the term describing the modulation of such oscillations, −𝛾𝜀 cos (𝛺𝑡 +𝜙)cos (𝜔𝑡)sin 𝜃, is the CC excitation. When the CC excitation is absent (𝜀= 0), we assume that the pendulum has a chaotic attractor for a given set of the remaining parameters. The chaotic behaviour of vertically excited pendula has been investigated for nearly forty years [24–29,31–35]. 2.1. Melnikov analysis Although the predictions from MM are both limited (only valid for motions based at points sufficiently near the separatrix) and approximate (the MM technique is a first-order perturbative method), they are of great interest due to the general scarcity of such analytical results in the theory of chaos. Since the technique has been described many times by different authors [36–38], we shall not discuss it in detail here, but analyse the results obtained from it. See Refs. [36–38] for more details about MM. Thus, it is also assumed that the complete system (1) satisfies the MM requirements, i.e., the excitation and dissipation terms are small-amplitude perturbations of the underlying integrable pendulum .. 𝜃+sin 𝜃= 0. After applying MM to Eq. (1), one straightforwardly obtains the Melnikov function (MF) 𝑀±(𝑡0)= −𝐷+𝐴sin (𝜔𝑡0)+𝐴′sin [(𝛺−𝜔)𝑡0+𝜙] +𝐴′′ sin [(𝛺+𝜔)𝑡0+𝜙], 𝐷≡8𝛿, 𝐴≡2𝜋𝛾𝜔2csch (𝜋𝜔∕2), 𝐴′≡𝜋𝜀𝛾 (𝛺−𝜔)2csch [𝜋(𝛺−𝜔)∕2], 𝐴′′ ≡𝜋𝜀𝛾 (𝛺+𝜔)2csch [𝜋(𝛺+𝜔)∕2],(2) where the positive (negative) sign refers to the top (bottom) homoclinic orbit of the conservative pendulum: 𝜃0(𝑡)= ± (4 arctan e𝑡−𝜋),. 𝜃0(𝑡)= ±2 sech 𝑡. To analyse the effect of the CC excitation, in the following we shall consider the normalized MF 𝑀± 𝑛(𝑡0)≡𝑀±(𝑡0)∕𝐷= −1 + 𝑅′sin [(𝛺−𝜔)𝑡0+𝜙] +𝑅sin (𝜔𝑡0)+𝑅′′ sin [(𝛺+𝜔)𝑡0+𝜙],(3) 𝑅≡𝐴∕𝐷, 𝑅′≡𝐴′∕𝐷, 𝑅′′ ≡𝐴′′∕𝐷. (4) The next two subsections study, for 𝑅 > 1(i.e., there exists homoclinic chaos in the absence of any CC excitation (𝜀= 0)), the appearance of simple zeros of 𝑀± 𝑛(𝑡0)for two subharmonic resonances: 𝛺=𝜔and 𝛺= 2𝜔. 2.1.1. Resonance 𝛺=𝜔 For the main resonance between the two excitations implicated, the normalized MF, Eq. (3), reduces to 𝑀± 𝑛(𝑡0)= −1 + 𝑅sin (𝜔𝑡0)+𝜀𝑅 sech (𝜋𝜔∕2)sin (2𝜔𝑡0+𝜙).(5) After expanding the factor sin (2𝜔𝑡0+𝜙)on the right-hand side of Eq. (5), one straightforwardly obtains 𝑀± 𝑛(𝑡0)= −1 − 𝜀𝑅 sech (𝜋𝜔∕2)sin 𝜙[2 sin2(𝜔𝑡0)− 1](6) +𝑅[1 + 2𝜀sech (𝜋𝜔∕2)cos 𝜙cos (𝜔𝑡0)] sin (𝜔𝑡0). Noting that min 𝑡0{1 + 𝜀𝑅 sech (𝜋𝜔∕2)sin 𝜙[2 sin2(𝜔𝑡0)− 1]} = 1 − 𝜀𝑅 sech (𝜋𝜔∕2)sin 𝜙, max 𝑡0[1+2𝜀sech (𝜋𝜔∕2)cos 𝜙cos (𝜔𝑡0)] = 1 + 2𝜀sech (𝜋𝜔∕2)cos 𝜙, and letting the CC excitation act on the pendulum such that 1+2𝜀sech (𝜋𝜔∕2)cos 𝜙⩽|𝜀sech (𝜋𝜔∕2)sin 𝜙− 1∕𝑅|,(7) one has that this relationship represents a sufficient condition for 𝑀± 𝑛(𝑡0)to be negative (or null) for all 𝑡0. The equals sign in Eq. (7) yields the approximate boundaries of the regions in the 𝜙−𝜀parameter plane in which homoclinic chaos is suppressed. After solving Eq. (7), one straightforwardly obtains the boundary functions 𝜀=(1 + 1∕𝑅)cosh (𝜋𝜔∕2) sin 𝜙− 2 cos 𝜙,(8) 𝜀=(1∕𝑅− 1)cosh (𝜋𝜔∕2) sin 𝜙+ 2 cos 𝜙.(9) One finds that 𝜙=𝜙sup,𝐼 𝑜𝑝𝑡,𝛺=𝜔≡2.68 and 𝜙=𝜙sup,𝐼𝐼 𝑜𝑝𝑡,𝛺=𝜔≡3.60 are the optimal suppressory initial phase differences for the boundaries (8) and (9), respectively, while the suppressory effectiveness of the CC excitation presents sensitivity to the specific initial chaotic state (i.e., to the value of the ratio 𝑅;Fig. 1 shows an illustrative example). Consider now the enhancing effect of the CC excitation. Noting that max 𝑡0{1 + 𝜀𝑅 sech (𝜋𝜔∕2)sin 𝜙[2 sin2(𝜔𝑡0)− 1]} = 1 + 𝜀𝑅 sech (𝜋𝜔∕2)sin 𝜙, min 𝑡0[1+2𝜀sech (𝜋𝜔∕2)cos 𝜙cos (𝜔𝑡0)] = 1 − 2𝜀sech (𝜋𝜔∕2)cos 𝜙 (cf. Eq. (6)), and letting the CC excitation act on the pendulum such that 1−2𝜀sech (𝜋𝜔∕2)cos 𝜙⩾|𝜀sech (𝜋𝜔∕2)sin 𝜙+ 1∕𝑅|,(10) then this relationship represents a sufficient condition for 𝑀± 𝑛(𝑡0)to change sign at some 𝑡0. Thus, enhancement of the initial homoclinic chaos is achieved by moving the pendulum from an effective homoclinic tangency condition (𝐷𝑒𝑓𝑓 =𝐴𝑒𝑓𝑓 )even more than in the initial situation with no CC excitation in the sense that 𝐷 < 𝐷𝑒𝑓 𝑓 ≡𝐷[𝜀sech (𝜋𝜔∕2)sin 𝜙+ 1] and 𝐴>𝐴𝑒𝑓𝑓 ≡𝐴[1−2𝜀sech (𝜋𝜔∕2)cos 𝜙] Physica D: Nonlinear Phenomena 457 (2024) 133950 3 R. Chacón et al. Fig. 1. Boundary functions given by Eqs. (8) ((b), (d)) and (9) ((a), (c)) encircling the regularization regions where homoclinic bifurcations are frustrated in the 𝜙−𝜀 parameter plane for the main resonance 𝛺=𝜔= 1.6and two values of the ratio 𝛿∕𝛾 (and hence of the ratio 𝑅, cf. Eqs. (2) and (3)): (a), (b) 𝑅= 1.1, and (c), (d) 𝑅= 2.852. The quantities plotted are dimensionless. over certain ranges of 𝜀and 𝜙. The equals sign in Eq. (10) yields the approximate boundaries of the regions in the 𝜙−𝜀parameter plane in which homoclinic chaos is enhanced. After solving Eq. (10), one straightforwardly obtains the boundary functions 𝜀=(1 + 1∕𝑅)cosh (𝜋𝜔∕2) 2 cos 𝜙− sin 𝜙,(11) 𝜀=(1 − 1∕𝑅)cosh (𝜋𝜔∕2) sin 𝜙+ 2 cos 𝜙.(12) One finds that 𝜙=𝜙enh,𝐼 𝑜𝑝𝑡,𝛺=𝜔≡5.81 and 𝜙=𝜙enh,𝐼𝐼 𝑜𝑝𝑡,𝛺=𝜔≡0.45 are the optimal enhancing initial phase differences for the boundaries (11) and (12), respectively, while the enhancing effectiveness of the CC excitation presents sensitivity to the specific initial chaotic state (i.e., to the value of the ratio 𝑅; see Fig. 2 for an illustrative example). Because of the perturbative nature of MM, the aforementioned values of the optimal suppressory and enhancing initial phase differences are approximate. It is worth noting that such optimal values present the property: |||𝜙sup,𝐼 𝑜𝑝𝑡,𝛺=𝜔−𝜙enh,𝐼 𝑜𝑝𝑡,𝛺=𝜔|||=|||𝜙sup,𝐼𝐼 𝑜𝑝𝑡,𝛺=𝜔−𝜙enh,𝐼𝐼 𝑜𝑝𝑡,𝛺=𝜔|||=𝜋. (13) 2.1.2. Resonance 𝛺= 2𝜔 In this case, the normalized MF, Eq. (3), can be recast in the form 𝑀± 𝑛(𝑡0)= −1 + 𝑅sin (𝜔𝑡0)+(1∕2)𝜀𝑅 sin (𝜔𝑡0+𝜙) +(9∕2)𝜀𝑅 1 + 2 cosh (𝜋𝜔)sin (3𝜔𝑡0+𝜙).(14) When 9[1 + 2 cosh (𝜋𝜔)]−1 ≪1, i.e., 𝜔 ≫ arccosh (4)∕𝜋≃ 0.656812, one can drop the third harmonic term in Eq. (14) and hence 𝑀± 𝑛(𝑡0)≲−1 + 𝑅√𝜀2∕4 + (1 + 𝜀cos 𝜙∕2)2.(15) If one now lets the CC excitation act on the pendulum such that 𝑅2(𝜀2∕4 + 1 + 𝜀2cos2𝜙∕4 + 𝜀cos 𝜙)− 1 ⩽0,(16) Fig. 2. Boundary functions given by Eqs. (11) ((b), (d)) and (12) ((a), (c)) encircling the regions where homoclinic chaos is enhanced in the 𝜙−𝜀parameter plane for the main resonance 𝛺=𝜔= 1.6and two values of the ratio 𝛿∕𝛾(and hence of the ratio 𝑅, cf. Eqs. (2) and (3): (a), (b) 𝑅= 1.1, and (c), (d) 𝑅= 2.852. The quantities plotted are dimensionless. this relationship represents a sufficient condition for 𝑀± 𝑛(𝑡0)to be negative (or null) for all 𝑡0. The equals sign in Eq. (16) yields the boundary of the region in the 𝜙−𝜀parameter plane in which homoclinic chaos is suppressed, 𝜀= − 2 cos 𝜙 1 + cos2𝜙±2√cos2𝜙−𝑅2+ 1 𝑅(1 + cos2𝜙),(17) with the constraint 𝑅 < √2, and where the sign +(−)before the square root corresponds to the upper (lower) branch of the boundary. One obtains that 𝜙=𝜙sup 𝑜𝑝𝑡,𝛺=2𝜔≡𝜋is the single optimal suppressory phase difference and that the effectiveness of the CC excitation noticeably depends on the specific initial chaotic state (value of the ratio 𝑅). Fig. 3 shows an example of the regions in the 𝜙−𝜀parameter plane in which homoclinic bifurcations are frustrated (Eqs. (16) and (17)). Also, one finds that the area of the regularization regions in the 𝜙−𝜀parameter plane for a fixed value of 𝑅, which provides a simple quantifier of the suppressory effectiveness of the CC excitation, is comparable for both resonances 𝛺={𝜔, 2𝜔}(compare the corresponding versions of Figs. 1 and 3). Consider now the enhancing effect of the CC excitation, also for the case 𝜔 ≫ arccosh (4)∕𝜋≃ 0.656812. Clearly, enhancement of the initial homoclinic chaos can be achieved by moving the pendulum from the homoclinic tangency condition (𝐷=𝐴)even more than in the initial situation with no CC excitation. This means that constraining −1 + 𝑅sin (𝜔𝑡0)to be in phase with (1∕2)𝜀𝑅 sin (𝜔𝑡0+𝜙)is a sufficient condition for the MF (14) to change sign at some 𝑡0. This condition provides the optimal enhancing initial phase difference, 𝜙=𝜙enh 𝑜𝑝𝑡,𝛺=2𝜔≡ 0, in the sense that the MF (14) presents its greatest maximum at 𝜙enh 𝑜𝑝𝑡,𝛺=2𝜔, i.e., one obtains the maximal gap from the homoclinic tangency condition. Also, the optimal suppressory and enhancing initial phase differences present the property |||𝜙sup 𝑜𝑝𝑡,𝛺=2𝜔−𝜙enh 𝑜𝑝𝑡,𝛺=2𝜔|||=𝜋, (18) as in the case of the main resonance. Physica D: Nonlinear Phenomena 457 (2024) 133950 4 R. Chacón et al. Fig. 3. Boundary functions (cf. Eq. (17)) encircling the regularization regions where homoclinic bifurcations are frustrated in the 𝜙−𝜀parameter plane for 𝜔= 1.6, 𝛺 = 2𝜔, and two values of the ratio 𝛿∕𝛾(and hence of the ratio 𝑅, cf. Eqs. (2) and (3)): (a) 𝑅= 1.1, (b) 𝑅= 1.35. The quantities plotted are dimensionless. 2.2. Numerical results We now contrast the analytical predictions derived from MM with the Lyapunov exponent (LE) calculations and bifurcation diagrams of the angular velocity for the parametrically excited damped pendulum (PEDP). Computer simulations of Eq. (1) showed that MM-based predictions supply meaningful qualitative information on the control scenario, even when the excitation amplitudes do not reasonably satisfy the MM requirements [36–38]. In this respect, it is worth recalling that even in the case of small amplitudes one cannot expect too good a quantitative agreement between these two types of approaches, because MM is a perturbative technique generally related to transient chaos, while LE supplies information solely about steady states. However, we shall show in this communication that MM-based predictions supply meaningful information about both suppression and enhancement of chaos for the present case of parametrically excited damped pendula under multiplicative control. To illustrate our findings, we shall use bifurcation diagrams which were constructed by means of a Poincaré map at the parameters indicated in the captions to the figures. Starting at a certain value of the control parameter of interest (for instance, 𝜀= 0), and taking the transient time as 5 × 102periods of the CI excitation after each of its increments (typically, we took 1.5 × 103values in the control parameter range), we sampled 103periods of the CI excitation by picking up the first . 𝜃values of every pulse cycle. For each value of the control parameter, we initialized on the last state found for the previous, slightly smaller, parameter value (i.e., “following the attractor”), while we typically chose the unstable equilibrium (𝜃=𝜋, 𝜃= 0) as the initial condition for its starting value. We computed LEs using a version of the algorithm introduced in [39], with integration typically up to 2 × 104drive cycles for a set of fixed parameters corresponding to an unavoidable situation (see Fig. 3 top in Ref. [19]) which is clearly outside the perturbative requirements (𝛿= 0.1, 𝛾 = 0.87, 𝜔 = 1.6). In the absence of any CC excitation (𝜀= 0), the PEDP exhibits a chaotic strange attractor characterized by a maximal LE 𝜆+(𝜀= 0)≃ 0.075 b/s. The maximal LE was calculated for each point on a uniform grid in the parameter plane 𝜙−𝜀. 2.2.1. Resonance 𝛺=𝜔 Fig. 4 shows the maximal LE distribution for the main resonance between the two excitations implicated. One sees that complete regularization (𝜆+(𝜀 > 0)<0)mainly appears in a region which corresponds to the theoretically predicted area but over the complete range of the modulation amplitude 𝜀(cf. Figs. 1(b) and 1(d)). For sufficiently small values of 𝜀, one finds enhancement of chaos over certain ranges of 𝜙which are roughly centred at the predicted values 𝜙=𝜙enh,𝐼 𝑜𝑝𝑡,𝛺=𝜔≡ 5.81 and 𝜙=𝜙enh,𝐼𝐼 𝑜𝑝𝑡,𝛺=𝜔≡0.45. Far from the perturbative regime, the Fig. 4. Maximal LE 𝜆+versus modulation amplitude 𝜀and phase difference 𝜙for 𝛿= 0.1, 𝛾 = 0.87, 𝜔 =𝛺= 1.6, and a grid of 75 ×75 points in the parameter plane 𝜙−𝜀. The quantities plotted are dimensionless. Fig. 5. Bifurcation diagrams of the angular velocity as a function of (a) the phase difference for 𝜀= 0.2and (b) the relative amplitude for 𝜙= 1. Inset in version (a): Enlargement of the window in which a period-doubling cascade occurs. The remaining parameters are the same as in Fig. 4. The quantities plotted are dimensionless. LE calculations indicate that maximum enhancement of chaos solely occurs over a narrow region roughly centred at 𝜙= 3𝜋∕2.Fig. 5 shows illustrative examples of the regularization routes as 𝜀and 𝜙are changed. Typically, the PEDP goes from the equilibrium 𝜃= . 𝜃= 0, which exists over a certain range starting at 𝜙= 0, to a strange chaotic attractor and, after an inverse period-doubling route, which is preceded by an inverse interior crisis, then to one of the two coexisting period-2 attractors, which exist over a wide range roughly centred at 𝜙=𝜙sup,𝐼 𝑜𝑝𝑡,𝛺=𝜔≡2.68, as the phase difference increases for a sufficiently large fixed value of 𝜀. Further increase of 𝜙leads to chaos via an infinite sequence of period-doubling bifurcations followed by an interior crisis (see the inset in Fig. 5(a)), with the chaotification route exhibiting approximate mirror symmetry with respect to the optimal suppressory value 𝜙sup,𝐼 𝑜𝑝𝑡,𝛺=𝜔≡2.68 (see Fig. 5(a)), and then to the destruction of this strange chaotic attractor via a boundary crisis, with the equilibrium 𝜃= . 𝜃= 0 being the single existing attractor from a certain value of the initial phase difference onwards (cf. Fig. 5(a)). Also, for a fixed phase difference 𝜙such that the PEDP presents a strange chaotic attractor at 𝜀= 0, this attractor is destroyed by a boundary crisis, again with the equilibrium 𝜃= . 𝜃= 0 being the single existing attractor from a certain value of the modulation amplitude onwards (cf. Fig. 5(b)). Physica D: Nonlinear Phenomena 457 (2024) 133950 5 R. Chacón et al. Fig. 6. Maximal LE 𝜆+versus phase difference 𝜙and modulation amplitude 𝜀for 𝛿= 0.1, 𝛾 = 0.87, 𝜔 = 1.6, 𝛺 = 2𝜔, and a grid of 75 ×75 points in the parameter plane 𝜙−𝜀. The quantities plotted are dimensionless. 2.2.2. Resonance 𝛺= 2𝜔 In this case, complete regularization (𝜆+(𝜀 > 0)<0)mainly appears in a region which approximately coincides with the theoretically predicted area (cf. Fig. 3), even for values of the modulation amplitude 𝜀clearly beyond the perturbative regime, as is shown in Fig. 6. One sees however that the symmetry of the regularization region is solely approximate from certain values of 𝜀onwards and that this area starts at much smaller 𝜀values than those predicted from MM (cf. Figs. 3 and 6). Fig. 7 shows illustrative examples of the regularization routes as 𝜀 and 𝜙are changed. Typically, the PEDP goes from a strange chaotic attractor at 𝜙= 0 to one of the two coexisting period-2 attractors at 𝜙=𝜙sup 𝑜𝑝𝑡,𝛺=2𝜔≡𝜋as the phase difference increases for a (not necessarily small) fixed value of the modulation amplitude (see Fig. 7(a)). The overall evolution of the initial chaotic state is characterized by the angular velocity undergoing an inverse period-doubling route as the phase difference is increased from 0, which is preceded by an inverse interior crisis. The regularization route presents approximate mirror symmetry with respect to the optimal suppressory value 𝜙sup 𝑜𝑝𝑡,𝛺=2𝜔≡𝜋. Also, even for a fixed phase difference relatively far from 𝜙sup 𝑜𝑝𝑡,𝛺=2𝜔≡𝜋, such as that considered in Fig. 7(b) (𝜙= 1.5), one finds a wide range of values of the modulation amplitude in which the pendulum’s dynamics is regularized. Note that the pendulum is chaotic for 𝜀= 0 while it is already regularized for 𝜀≈ 5 × 10−4 (see Fig. 7(b)). 2.2.3. Other resonances As suggested by comparison of the above findings for the resonances 𝛺=𝜔and 𝛺= 2𝜔, one could guess that multiplicative control of the PEDP would allow an even wider range of effective resonances. We numerically found that this is indeed the case. Fig. 8 shows the maximal LE distribution for the exemplary case of the resonance 𝛺=𝜔∕2. In this case, complete regularization (𝜆+(𝜀 > 0)<0)mainly appears over certain ranges of 𝜙which are roughly centred at the values 𝜙= 𝜙sup,𝐼 𝑜𝑝𝑡,𝛺=𝜔∕2 ≡𝜋∕2 and 𝜙=𝜙sup,𝐼𝐼 𝑜𝑝𝑡,𝛺=𝜔∕2 ≡3𝜋∕2, respectively, for a wide range of 𝜀values, while enhancement of chaos occurs over certain ranges of 𝜙which are roughly centred at the values 𝜙=𝜙𝑒𝑛ℎ,𝐼 𝑜𝑝𝑡,𝛺=𝜔∕2 ≡0 and 𝜙=𝜙𝑒𝑛ℎ,𝐼𝐼 𝑜𝑝𝑡,𝛺=𝜔∕2 ≡𝜋for sufficiently large values of 𝜀. Remarkably, such optimal values present the property ||||𝜙sup,𝐼 𝑜𝑝𝑡,𝛺=𝜔 2 −𝜙enh,𝐼 𝑜𝑝𝑡,𝛺=𝜔 2||||=||||𝜙sup,𝐼𝐼 𝑜𝑝𝑡,𝛺=𝜔 2 −𝜙enh,𝐼𝐼 𝑜𝑝𝑡,𝛺=𝜔 2||||=𝜋 2.(19) Fig. 9 shows that the regularization routes as the parameters 𝜀and 𝜙 are independently changed are similar to those found for the resonances Fig. 7. Bifurcation diagrams of the angular velocity as a function of (a) the initial phase difference for 𝜀= 0.2and (b) the modulation amplitude for 𝜙= 1.5. Insets: Enlargement of the windows in which inverse and direct period-doubling cascades occur in versions (a) and (b), respectively. The remaining parameters are the same as in Fig. 6. The quantities plotted are dimensionless. Fig. 8. Maximal LE 𝜆+versus phase difference 𝜙and modulation amplitude 𝜀for 𝛿= 0.1, 𝛾 = 0.87, 𝜔 = 1.6, 𝛺 =𝜔∕2, and a grid of 75 ×75 points in the parameter plane 𝜙−𝜀. The quantities plotted are dimensionless. Fig. 9. Bifurcation diagrams of the angular velocity as a function of (a) the initial phase difference for 𝜀= 0.2and (b) the modulation amplitude for 𝜙= 1. The remaining parameters are the same as in Fig. 8. The quantities plotted are dimensionless. 𝛺={𝜔, 2𝜔}. But the regularization (enhancement) ranges are narrower (wider) than those corresponding to resonances 𝛺={𝜔, 2𝜔}(compare the corresponding versions of Figs. 5 and 9). Fig. 9(b) shows that the pendulum is chaotic for 𝜀= 0 while it is already regularized for 𝜀≈ 5 × 10−4through an inverse period doubling cascade. Physica D: Nonlinear Phenomena 457 (2024) 133950 6 R. Chacón et al. 3. Controlling chaos in starlike networks Previous studies have shown that chaos in coupled arrays of damped, periodically driven, nonlinear oscillators can be suppressed by random shortcuts [40], parametric disorder [41], global disordered driving forces [42], impurities [43], and localized resonant forces [18]. Also, in the context of networks of nonchaotic (overdamped) systems, such as periodically driven bistable oscillators, diverse amplification phenomena have been studied in the presence [44] and in the absence of noise [45]. Recently, the problem of how network topology affects a system’s controllability has been considered using a canonical linear time-invariant dynamics [46]. We consider in the following the application of the above control scenario to a topology consisting of a starlike network (one hub and 𝑁− 1 peripheral leaves) where each node is occupied by a chaotic pendulum and where the CC excitation is solely applied to a number, 𝑀, of pendula: .. 𝜃𝑖+ sin 𝜃𝑖= −𝛾[1 + 𝜎𝐻𝐹𝑐(𝑡)]cos (𝜔𝑡)sin 𝜃𝑖 −𝛿 . 𝜃𝑖+𝜆sin (𝜃𝐻−𝜃𝑖), .. 𝜃𝐻+ sin 𝜃𝐻= −𝛿 . 𝜃𝐻−𝛾[1 + 𝜎𝑖𝐹𝑐(𝑡)]cos (𝜔𝑡)sin 𝜃𝐻 +𝜆 𝑁−1 ∑ 𝑖=1 sin (𝜃𝑖−𝜃𝐻),(20) 𝑖= 1,…, 𝑁 − 1. These equations describe the dynamics of a highly connected node (or hub), 𝜃𝐻, and 𝑁− 1 linked pendula (or leaves), 𝜃𝑖, with 𝐹𝑐(𝑡)≡𝜀cos (𝛺𝑡 +𝜑)being the (local) CC excitation, while 𝜎𝐻 (𝜎𝑖) is equal to 1 when the CC excitation acts on the hub (leaf 𝑖) and 0 otherwise, while 𝜆is the coupling constant. Thus, with a fixed size of the network, 𝑁, fixed values of the parameters 𝛿, 𝛾, and 𝜔corresponding to a given initial chaotic state of the uncoupled pendula, and fixed suitable (suppressory) values of 𝜀, 𝜙, and 𝛺(cf. Sec. II), we study the network regularization in the regime of both weak and strong coupling (0< 𝜆 ⩽2)as a function of two independent parameters, 𝑀and 𝜆. The dynamic equation of these networks [Eq. (20)] was numerically integrated using a fourth-order Runge–Kutta algorithm, while the initial conditions were chosen randomly and independently for each node of the networks. To visualize the global spatiotemporal dynamics of the networks, we calculate the average velocity 𝜎(𝑗𝑇 )≡1 𝑁 𝑁 ∑ 𝑛=1 . 𝜃𝑛(𝑗𝑇 ),(21) where 𝑗is an integer multiple of the excitation period 𝑇≡2𝜋∕𝜔, while the degree of synchronization is characterized by the correlation function [14,47] 𝐶≡2 𝑁(𝑁− 1) ∑ (𝑖𝑙)⟨cos (𝜃𝑖−𝜃𝑙)⟩𝑡,(22) with the summation being over all pairs of pendula, and where ⟨⋅⟩𝑡indicates time averaging over a predefined (sufficiently long) observation window. Note that 𝐶is 1for the perfectly synchronized state while desynchronization increases as 𝐶decreases from 1. Let us consider the effect of the multiplicative control on a single peripheral pendulum 𝜃𝑗(𝑀= 1)while the remaining pendula, including the hub, are free from control. Note that, a priori, this could be the least favourable case in terms of completely regularizing the whole network. Numerical simulations indicate, however, that regularization to periodic states is possible over certain coupling intervals even for relatively small 𝜀values, although the scenario is strongly dependent on the initial conditions. For the sake of clarity, consider first the case where all pendula are uncoupled (𝜆= 0)and present a perfectly synchronized chaotic state (scenario I). In this situation, one can apply the control excitation to an isolated (peripheral when coupled) pendulum regularizing its behaviour, then increase the coupling parameter to a small positive value 𝜆 ≳ 0and determine the new stationary state of such a peripheral Fig. 10. Bifurcation diagrams of (a) the average angular velocity, (b) the correlation function, and (c) the maximal LE 𝜆+as functions of the coupling constant 𝜆, while a single peripherical pendulum is subjected to multiplicative control (node 1). (d) Angular velocities versus time for 𝜆= 0.22 and (e) 𝜆= 0.8.𝜀= 0.2,𝜙= 1,𝑁= 10, while the hub corresponds to node 10. The remaining parameters are the same as in Fig. 4. The quantities plotted are dimensionless and all pendula are synchronized at 𝑡= 0. pendulum (and of all of the remaining pendula of the starlike network). After using this new state as initial condition, one further increases the coupling parameter and follows the attractors of all the pendula as 𝜆is increased. In general, we found that complete synchronization of the pendula is not possible (see Figs. 10(b), 10(d), and 10(e)), while, even for relatively small values of the coupling parameter, the single peripheral pendulum can regularize the complete network, as shown in Figs. 10(a), 10(c), and 10(d). This latter figure corresponds to a set of parameters for which an isolated pendulum is chaotic in the absence of any control excitation, while its behaviour is periodic for 𝜀= 0.2,𝜙= 1 (cf. Figs. 4 and 5). Let us consider now a different case where all pendula are uncoupled (𝜆= 0)and present an asynchronous chaotic state due to the initial conditions of all pendula being independently and randomly chosen (scenario II). After applying the same coupling procedure as in scenario I, we typically found that the network remains asynchronous and chaotic even for relatively large values of the coupling constant, as is shown in Fig. 11(b), 11(c), and 11(d). Remarkably, there exists a Physica D: Nonlinear Phenomena 457 (2024) 133950 7 R. Chacón et al. Fig. 11. Bifurcation diagrams of (a) the average angular velocity, (b) the correlation function, and (c) the maximal LE 𝜆+as functions of the coupling constant 𝜆while a single peripherical pendulum is subjected to multiplicative control (node 1). (d), (e) Angular velocities versus time for 𝜀= 0.5,𝜙= 0.68,𝑁= 10, the hub being node 10, and two values of the coupling constant: (d) 𝜆= 0.5and (e) 𝜆= 1.3. The remaining parameters are the same as in Fig. 4. The quantities plotted are dimensionless and all pendula are desynchronized at 𝑡= 0. critical coupling value, close to 𝜆= 1, for which the network reaches a perfectly synchronized chaotic state, as in the instance shown in Fig. 11(e), while for an even higher value of 𝜆all pendula regularize, as is shown in Fig. 11(a) and 11(c). In general, our numerical results confirmed that the MM-based analytical predictions for an isolated pendulum remain useful for the starlike network over significant ranges of the coupling constant. We also found that desynchronization states appear due to cluster synchronization of different sets of pendula over almost the entire range of the coupling constant (cf. Figs. 10 and 11) when the size of the network is large enough. Next, it is interesting to study the cumulative effect of applying the multiplicative control on an increasing number of peripheral pendula (𝑀 > 1), while the hub remains free from control, in the weak coupling regime where synchronization phenomena do not yet dominate the networks’ dynamics. We assume an initial situation where regularization is not possible for any value of the modulation amplitude 𝜀when the control is applied to a single peripheral pendulum (𝑀= 1)and all the pendula are Fig. 12. (a),(b), (c): Bifurcation diagrams of the average angular velocity 𝜎as a function of the modulation amplitude 𝜀corresponding to three cases of multiplicative control: (a) 𝑀= 1, (b) 𝑀= 6 and (c) 𝑀= 9. (d) Correlation function 𝐶as a function of the modulation amplitude 𝜀for 𝑀= 1 (solid blue line), 𝑀= 6 (dashed red line) and 𝑀= 9 (dotted black line). Fixed parameters: 𝑁= 10,𝜆= 0.1, 𝜙 = 1, 𝛿= 0.1, 𝛾 = 0.87, 𝜔 =𝛺= 1.6. The quantities plotted are dimensionless and all pendula are initially synchronized. synchronized. Following a similar procedure to that described above for scenario I, by increasing 𝑀from unity, one typically obtains regularization of the whole network for sufficiently large 𝜀values on the one hand, together with an improvement of the synchronization of the regularized dynamics on the other, as in the instance shown in Fig. 12 for 𝜆= 0.1 (weak coupling regime). By again increasing 𝑀from unity but with all pendula initially chaotic and asynchronous, we found a scenario similar to scenario II described above for a single peripheral pendulum. Specifically, synchronization and regularization occur from a certain value of the coupling constant, 𝜆=𝜆(𝑀), which decreases with the number 𝑀of nodes subjected to multiplicative control, as shown in Fig. 13, in which the parameters correspond to a situation where an isolated pendulum subjected to multiplicative control presents regular behaviour. Finally, we studied the role played by the degree of connectivity in the chaos-control scenario by applying the multiplicative control to the central pendulum. In the case of a single control, one finds that controlling the most highly connected pendulum is by far the most effective control procedure, as shown in Fig. 14. Remarkably, solely applying the multiplicative control to the hub (𝑀= 1)is a much better choice than controlling even several peripheral pendula but not the hub (compare the respective versions of Fig. 13(a)–(c) with those of Fig. 14). The reason for this relatively good effectiveness stems from two facts. First, solely controlling the hub does not significantly break the synchronization of the whole network when 𝑁is sufficiently large (compare Figs. 10(b) and 14(a)). Second, its maximum degree of connectivity allows the hub to directly influence all the remaining (peripheral) pendula –in the sense of taming their chaotic dynamics– due to its behaving as a direct ordering source for all of them for sufficiently large 𝜀values (compare Figs. 11(a) and 14(b)). Once the hub has been subjected to multiplicative control, one could expect a priori that additionally controlling other (peripheral) pendula should improve the network’s regularization. As expected, when all pendula are subjected to the same multiplicative control, the network’s synchronization becomes perfect and the regularization route as the modulation amplitude is varied coincides with that of an isolated pendulum subjected to the same remaining parameters, involving typically several consecutive crises followed by an inverse period doubling to finally reach the equilibrium (𝜃= 0, . 𝜃= 0)when the common amplitude modulation is sufficiently large. Similar results have been obtained for other resonances, as shown in Figs 15 and 16 for the resonances 𝛺= 2𝜔and 𝛺=𝜔∕2, respectively. In particular, one sees that the higher effectiveness of controlling only the hub instead of a single peripheral pendulum, in the sense of obtaining regularization of the starlike network for smaller values of the coupling, holds beyond the main resonance case (𝛺=𝜔). Physica D: Nonlinear Phenomena 457 (2024) 133950 8 R. Chacón et al. Fig. 13. Bifurcation diagrams of the average velocity 𝜎and correlation function 𝐶as functions of the coupling parameter 𝜆for 𝜀= 0.5,𝜙= 0.68,𝑁= 10, the hub being the node 10. (a) 𝑀= 1. (b) 𝑀= 2. (c) 𝑀= 3. (d) 𝑀= 4. (e) 𝑀= 5. (f) 𝑀= 6. (g) 𝑀= 7. (h) 𝑀= 8. (i) 𝑀= 9. The remaining parameters are the same as in Fig. 4. The quantities plotted are dimensionless and all pendula are initially desynchronized. 4. Concluding remarks We have shown theoretically and numerically that the application of suitable amplitude modulations is a reliable procedure to control (suppress and enhance) the chaotic behaviour of both isolated parametrically excited damped pendula and starlike networks of them subjected to sinusoidal coupling. Effective regularization of the entire network is typically obtained under localized control on an increasing number of pendula when the coupling is large enough. Remarkably, this reliable regularization has been demonstrated for several resonances between the two excitations implicated, being a genuine feature of the multiplicative control in contrast with additive control in which the main resonance is by far the most effective. We have shown how the effectiveness of this multiplicative control, when the amplitude modulation is applied to a single node, strongly depends upon the degree of the target node: applying the multiplicative control to the highest-degree node is by far the best suppressory strategy, while Fig. 14. Bifurcation diagrams of the average velocity 𝜎and correlation function 𝐶as functions of the coupling parameter 𝜆for 𝑁= 10. (a) All pendula are synchronized at 𝑡= 0,𝜀= 0.2,𝜙= 1. (b) All pendula are asynchronous at 𝑡= 0,𝜀= 0.5, 𝜙= 0.68. The hub (node 10) is the single pendulum subjected to multiplicative control, while the remaining parameters are the same as in Fig. 4. The quantities plotted are dimensionless. Fig. 15. Bifurcation diagrams of the average velocity 𝜎and correlation function 𝐶as functions of the coupling parameter 𝜆, while all pendula are asynchronous at 𝑡= 0. Multiplicative control applied to (a) a single peripherical node and (b) the hub for 𝛺= 2𝜔= 1.6,𝜀= 0.2,𝜙= 1.5, and the same remaining parameters as in Fig. 6. The quantities plotted are dimensionless. Fig. 16. Bifurcation diagrams of the average velocity 𝜎and correlation function 𝐶 as functions of the coupling parameter 𝜆, while all pendula are asynchronous at 𝑡= 0. Multiplicative control at (a) a single peripherical node and (b) the hub for 𝛺=𝜔∕2 = 0.8,𝜖= 0.2,𝜙= 1.0, and the same remaining parameters as in Fig. 9. The quantities plotted are dimensionless. Physica D: Nonlinear Phenomena 457 (2024) 133950 9 R. Chacón et al. applying it to low-degree nodes is the poorest choice. In the case of applying multiplicative control to several low-degree nodes, we have found a systematic improvement of the network’s regularization and synchronization. Finally, we hope that our results may be useful to optimally control chaos in scale-free networks of damped-driven oscillators since a highly connected node in such a network can be thought of as a hub of a locally starlike part of the network, with a degree of connectivity that belongs to the complete network’s degree-of-connectivity distribution. Other problems that deserve to be explored are, on the one hand, the robustness of the present multiplicative control against changes of the amplitude-modulation waveform and, on the other, the presence of time-varying coupling in temporal networks [48]. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgements Financial support from the Ministerio de Ciencia, Innovación y Universidades (MICIU, Spain) through Project No. PID2019-108508GB- I00/AEI/10.13039/501100011033 cofinanced by FEDER funds (R.C., A.M.G.-H., F.P.) and from the Junta de Extremadura (JEx, Spain) through Project No. GR21012 cofinanced by FEDER funds (R.C.) is gratefully acknowledged. References [1] G. Chen, X. Dong, From Chaos To Order World Scientific, Singapore, 1998. [2] S. Boccaletti, et al., Phys. Rep. 329 (2000) 103. [3] R. Chacón, Control of Homoclinic Chaos By Weak Periodic Perturbations, World Scientific, London, 2005. [4] K. Pyragas, Phys. Lett. A 170 (1992) 421. [5] W.L. Ditto, S.N. Rauseo, M.L. Spano, Phys. Rev. Lett. 65 (1990) 3211. [6] A. Azevedo, S.M. Rezende, Phys. Rev. Lett. 66 (1991) 1342. [7] E.R. Hunt, Phys. Rev. Lett. 67 (1991) 1953. [8] K. Murali, S. Sinha, Phys. Rev. E 68 (2003) 016210. [9] S. Zambrano, et al., Chaos 16 (2006) 013111. [10] P.J. Martínez, S. Euzzor, J.A.C. Gallas, R. Meucci, R. Chacón, Sci. Rep. 7 (2017) 17988. [11] A. Kandangath, S. Krishnamoorthy, Y.-C. Lai, J.A. Gaudet, IEEE Trans. Circuits Syst. 54 (2007) 1109. [12] Y. Braiman, I. Goldhirsch, Phys. Rev. Lett. 66 (1991) 2545. [13] J.M. Seoane, et al., Phys. Rev. E 78 (2008) 016205. [14] R. Chacón, F. Palmero, J. Cuevas-Maraver, Phys. Rev. E 93 (2016) 062210. [15] I.B. Schwartz, I. Triandaf, R. Meucci, T.W. Carr, Phys. Rev. E 66 (2002) 026213. [16] R. Meucci, et al., Physica D 189 (2004) 70. [17] A. Farshidianfar, A. Saghafi, Phys. Lett. A 378 (2014) 3457. [18] P.J. Martínez, R. Chacón, Phys. Rev. Lett. 93 (2004) 237006, 96, 059903(E) (2006). [19] R. Chacón, A. Martínez García-Hoz, J.A. Martínez, Phys. Rev. E 95 (2017) 052219. [20] V. Ravichandran, V. Chinnathambi, S. Rajasekar, Physica A 376 (2007) 223. [21] M. Siewe Siewe, C. Tchawoua, S. Rajasekar, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 21 (2011) 1583. [22] C.H. Miwadinou, A.V. Monwanou, L.A. Hinvi, J.B. Chabi Orou, Chaos Solit. Fract. 113 (2018) 89. [23] R. Chacón, A. Martínez García-Hoz, J.J. Miralles, P.J. Martínez, Phys. Lett. A 378 (2014) 1104. [24] S. Lenci, G. Rega, Physica D 240 (2011) 814. [25] R.W. Leven, B.P. Koch, Phys. Lett. A 86 (1981) 71. [26] M. Bartuccelli, P.L. Christiansen, N.F. Pedersen, M.P. Soerensen, Phys. Rev. B 33 (1986) 4686. [27] J. Yang, Z. Jing, Chaos Solit. Fract. 42 (2009) 1214. [28] X. Chen, Z. Jing, X. Fu, Nonlinear Dynam. 78 (2014) 317. [29] L. Zhou, F. Chen, Internat. J. Bifur. Chaos Appl. Sci. Engrg. 30 (2020) 2050166. [30] R. Milo, et al., Science 298 (2002) 824. [31] B.P. Koch, R.W. Leven, Physica D 16 (1985) 1. [32] J.A. Blackburn, N. Gronbeck-Jensen, H.J.T. Smith, Phys. Rev. Lett. 74 (1995) 908. [33] E. Butikov, Amer. J. Phys. 69 (2001) 755. [34] U. Lepik, H. Hein, J. Sound Vibr. 288 (2005) 275. [35] J. Sieber, et al., Phys, Rev. Lett. 100 (2008) 244101. [36] V.K. Melnikov, Trans. Moscow Math. Soc. 12 (1963) 1. [37] J. Guckenheimer, P.J. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, Berlin, 1983. [38] A.J. Lichtenberg, M.A. Lieberman, Regular and Stochastic Motion, Springer, New York, 1983. [39] G. Benettin, L. Galgani, J.M. Strelcyn, Phys. Rev. A 14 (1976) 2338; I. Shimada, T. Nagasama, Prog. Theor. Phys. 61 (1979) 1605. [40] F. Qi, Z. Hou, H. Xin, Phys. Rev. Lett. 91 (2003) 064102. [41] S.H. Strogatz, Nature 378 (1995) 444. [42] R. Chacón, P.J. Martínez, Phys. Rev. Lett. 98 (2007) 224102. [43] M. Weiss, T. Kottos, T. Geisel, Phys. Rev. E 63 (2001) 056211. [44] A.S. Pikovsky, A. Zaikin, M.A. de la Casa, Phys. Rev. Lett. 88 (2002) 050601. [45] J.A. Acebrón, S. Lozano, A. Arenas, Phys. Rev. Lett. 99 (2007) 128701; J. Zhou, Y. Zhou, Z. Liu, Phys. Rev. E 83 (2011) 046107. [46] Y.-Y. Liu, J.-J. Slotine, A.-L. Barabási, Nature 473 (2011) 167. [47] R. Chacón, A. Martínez García-Hoz, P.J. Martínez, J.A. Martínez, submitted for publication. [48] Y. Zhang, S.H. Strogatz, Nature Commun. 12 (2021) 3273.