# Desynchronization and pattern formation in a noisy feedforward oscillators network

###### Abstract

We consider a one-dimensional directional array of diffusively coupled oscillators. They are perturbed by the injection of a small additive noise, typically orders of magnitude smaller than the oscillation amplitude, and the system is studied in a region of the parameters that would yield deterministic synchronization. Non normal directed couplings seed a coherent amplification of the perturbation: this latter manifests as a modulation, transversal to the limit cycle, which gains in potency node after node. If the lattice extends long enough, the initial synchrony gets eventually lost and the system moves toward a non trivial attractor, which can be analytically characterized as an asymptotic splay state. The noise assisted instability, ultimately vehiculated and amplified by the non normal nature of the imposed couplings, eventually destabilizes also this second attractor. This phenomenon yields spatiotemporal patterns, which cannot be anticipated by a conventional linear stability analysis.

## 1 Introduction

Understanding the origin and functional significance of self-organized patterns of activity, is a challenging question of broad applied and fundamental importance [1, 2]. In many realms of investigation, the system under inspection is composed by individual excitable units, which execute periodic oscillations [3]. Often, coupling together an ensemble made of identical oscillators can eventually yield a fully synchronized solution [4]. This amounts to operate the system in unison, the oscillations displayed on different sites of the collection being perfectly coordinated, with no phase delay. For many applications of interest, as e.g. the study of collective oscillations in neuroscience, distinct deterministic oscillators occupy the nodes of a heterogenous network, which defines the embedding structural support [5, 6]. Diffusive couplings between adjacent mesoscopic units are customarily assumed, a paradigmatic choice which proves adequate in many cases [7], from modeling the electrical synapses to problems related to the energy management in power plants. Instabilities, however, may be triggered by the punctual injection of a heterogeneous perturbation [8], a tiny source of stochastic disturbance which, under specific conditions, amplifies and eventually breaks the oscillators’ synchrony [9]. The instabilities instigated by random fluctuations are often patterns precursors [10, 11]. The imposed perturbation materializes in fact in patchy motifs of the concentration amount, characterized by a vast gallery of shapes and geometries. An archetypal model of self-sustained oscillations is the celebrated Complex Ginzburg-Landau equation (CGLE), often evoked as a pillar of non linear phenomena, from superconductivity to superfluidity and Bose-Einstein condensation, via strings in field theory and neuroscience [12]. The CGLE, defined on ordinary or graph-like supports, admits a time-dependent uniform synchronized solution, of the limit cycle type. Deviations from a periodic waveform, sustained by nonlinearities, yield a prototypical modulational instability characterized by spectral-sidebands and the breakup of the waveform into a train of pulses. This is the so called Benjamin-Feir (BF) instability, named after the researchers who first identified the phenomenon working with periodic surface gravity waves (Stokes waves) on deep water [13]. Typically the condition for the onset of the deterministic instability can be straightforwardly worked out by means of a traditional linear stability analysis, which constraints the reaction parameters involved in the formulation of the problem [1, 14, 15, 16].

Starting from these premises we are here interested in studying the stochastic analog of the BF instability in an open feed-forward topology. In the framework that we shall set to explore, the complex state variable of the CGLE is disturbed by a small exogenous perturbation, which configures as an additive white noise, possibly orders of magnitude smaller than the unperturbed oscillation amplitude. More specifically, we are interested in assessing the role played by the injected stochastic drive, when the system is operated in a parameters region for which the synchronous limit cycle proves stable under deterministic evolution. As we shall argue, synchronous solutions, deemed deterministically stable, can turn unstable by agitating the system with an arbitrarily small perturbation. In the case of the CGLE, here assumed as a reference model, we will unfold, and thoroughly characterize, a generalized class of convective instabilities reminiscent of the BF one. To achieve the sought effect we shall accommodate for non-normal [17], diffusive couplings between individual oscillators. In a recent series of papers [18, 19], we showed that giant stochastic oscillations, with tunable frequencies, can be obtained, by replicating a minimal model for quasi-cycle along a directed chain of coupled oscillators. Here, the directed link between adjacent oscillators will fuel a self-consistent amplification of the stochastic disturbance, always yielding – for a sufficiently long chain – a loss of synchronicity. Taken all together, our findings constitute a practical example of convective instability [20], and point to the subtle interplay between noise and topology, capable of changing qualitatively the system dynamics. This brings into evidence possible failure of a purely deterministic approaches to real life problems.

The paper is organised as follows: in the next section we will introduce the model to be probed. In particular, we will discuss its synchronized and splay states. We shall then turn to analyze the effect of stochasticity, with reference to the amplification mechanism, as alluded to above. Stochastic non normal patterns are consequently reported to occur, notwithstanding the stability of the homogeneous solution under traditional linear analysis. Finally, we will sum up and draw our conclusions.

## 2 Deterministic Ginzburg-Landau oscillators: synchronized and splay states

Our model consists of \Omega diffusively and unidirectionally coupled Ginzburg-Landau oscillators. Each oscillator is described by the complex variable W_{j} (1\leq j\leq\Omega). The oscillators in this directionally coupled chain (see Fig.1) obey the following ordinary differential equations

\frac{\mathrm{d}W_{1}}{\mathrm{d}t}=W_{1}-(1+ic_{2})|W_{1}|^{2}W_{1} | (1a) | ||

and, for j>1 | |||

\frac{\mathrm{d}W_{j}}{\mathrm{d}t}\!\!=\!\!W_{j}\!-\!(1\!+ic_{2})|W_{j}|^{2}W% _{j}+(1+ic_{1})K(W_{\!j-1}\!-\!W_{j}) | (1b) |

where c_{1}, c_{2} are real parameters and K denotes the coupling strength. It is obvious that changing the sign of K and at the same time inverting the boundary conditions is equivalent to reversing the information flow along the chain: therefore in the rest of this paper K is assumed to be positive. The system is also symmetric under the following transformation: W_{j}\rightarrow W_{j}^{*}, (c_{1},c_{2})\rightarrow-(c_{1},c_{2}) which allows us to restrict our focus on half of the (c_{1},c_{2}) parameter plane. Two types of solution are of interest, the synchronized and the splay ones. The synchronized state (usually denoted as homogeneous state, in the vast literature of spatially coupled oscillators) corresponds to the solution

W_{j}=\exp(-ic_{2}t)\quad,\quad j=1,\cdots,\Omega. | (2) |

By direct inspection of Eq. (1) one can check that any dependence on the spatial coupling K and on the the parameter c_{1} disappears, and that this solution exists for any value of c_{2}.

The splay states constitute a family of uniformly rotating solutions with finite constant-in-time phase differences between consecutive nodes. Inserting the general polar form W_{j}=\rho_{j}\exp(-i\theta_{j}) into equation (1), and assuming that \dot{\rho_{j}}=0 and \theta_{j}=-c_{2}t+\sum^{j}_{k=1}\phi_{k} (where the \phi_{k}=\theta_{k}-\theta_{k-1} are constant-in-time phase differences) leads to the recurrence relations

\rho_{j}=\sqrt{\left(1+K\left[\frac{\rho_{j-1}}{\rho_{j}}f(\phi_{j})-1\right]% \right)} | (3a) | ||

0=c_{2}(1-\rho_{j}^{2})+K\left[\frac{\rho_{j-1}}{\rho_{j}}g(\phi_{j})-c_{1}\right] | (3b) |

where

f(\phi_{j})=\cos\phi_{j}+c_{1}\sin\phi_{j} | (4a) | ||

g(\phi_{j})=c_{1}\cos\phi_{j}-\sin\phi_{j}. | (4b) |

The recurrence relations (3) need to be completed by a suitable initial condition, where the state variables in the first node of the chain are fixed to some values: we adopt the condition \rho_{1}=1 and \phi_{1}=0, that would be compatible with the synchronous solution (2). We avoid reporting explicit calculations, but it can be shown that the recurrence relations beyond the homogeneous solution, also admits a second solution with non zero \phi_{j}, which spatially converges to the asymptotic splay state

\rho_{\infty}=\sqrt{1+K\left(f(\phi_{\infty})-1\right)}\\ | (5a) | ||

\phi_{\infty}=2Atan\left[\frac{1+c_{1}c_{2}}{c_{2}-c_{1}}\right] | (5b) |

where the special case \phi_{\infty}=\pm\pi occurs in the limit c_{2}\rightarrow c_{1}. In practice, one finds that the asymptotic splay state is rapidly approached along the chain (see Fig.(2))

The rate of convergence depends on the paramters K, c_{1} and c_{2}, but here, for the sake of space, we do not report any detailed investigation on this point.

It is important to point out that the existence condition for the splay state is that \rho_{\infty} is real, i.e. that the argument of the square root in Eq. (5a) is non-negative. As an example, in Fig. 3 we show the region in the (c_{1},c_{2})-plane where the splay state exists for K=4: the colour code corresponds to different positive values of \rho_{\infty}, while the black region indicates where the splay state does not exist.

As a final remark, we want to point out that there exist an entire family of solutions asymptotically approaching along the chain the splay state (see Eqs. (3)). In these solutions the synchronous state extends to an arbitrary large initial portion of the chain, namely \rho_{j}=1 and \phi_{j}=0 for j=1,\cdots,{\bar{j}}. For j>\bar{j} phase differences become finite and the solution converges to the asymptotic splay state (5) for large j. As we shall discuss later, the existence of this entire family of splay states impacts on the way noise destabilizes the homogeneous synchronized state determining a typical spatio-temporal pattern organization for the stochastic system.

### 2.1 Stability of synchronized and splay states

In order to investigate the stability of the synchronous and of the splay states we can perform a standard linear stability analysis. We first introduce small perturbations \delta\rho_{j},\delta\theta_{j} of the limit cycles, W_{j}=(\rho_{j}+\delta\rho_{j})\exp(i(\theta_{j}+\delta\theta_{j})) for 1\leq j\leq\Omega. Linearizing and retaining the first order in the pertubations leads to an equation that can be put in the general matrix form \delta\dot{\bf v}={\bf J}({\bf\rho},{\bf\theta})\delta{\bf v} (we adopt the shorthand notation {\bf(\rho,\phi)}=(\rho_{1},\phi_{1},\rho_{2},\phi_{2},\cdots,\rho_{\Omega},% \phi_{\Omega})), where \delta{\bf v}=\delta(\rho,\phi) is the vector of perturbations and {\bf J} is the Jacobian matrix associated to dynamics (1). Due to the unidirectional nature of the coupling K, {\bf J} exhibits a lower tridiagonal block structure. Hence, to assess the stability of any state it is enough to compute the eigenvalues \lambda_{\rho_{j}} and \lambda_{\theta_{j}} for 1\leq j\leq\Omega of the diagonal 2\times 2 blocks

{\bf A}_{1}=\left(\begin{array}[]{lr}-2&\;\;0\\ &\\ -2c_{2}&\;\;0\end{array}\right) | (6a) | ||

and | |||

{\bf A}_{j}\!\!=\!\!\left(\begin{array}[]{ll}\!(1-3\rho_{j}^{2}-K)&\;\;K\rho_{% j-1}g(\phi_{j})\\ &\\ \!-\!\left(2c_{2}\rho_{j}+K\frac{\rho_{j-1}}{\rho_{j}}g(\phi_{j})\right)&\;\;-% K\frac{\rho_{j-1}}{\rho_{j}}f(\phi_{j})\end{array}\!\!\right) | (6b) |

for 2\leq j\leq\Omega.

The eigenvalues of the first block A_{1} are \lambda_{\rho_{1}}\!\!=\!\!-2 and \lambda_{\theta_{1}}\!\!=\!\!0, the latter reflecting marginal stability towards global phase rotations. A given limit cycle solution is stable only if the all other block complex eigenvalues have a negative real part, i.e. \Re(\lambda_{\rho_{j}})<0 and \Re(\lambda_{\theta_{j}})<0 for 2\leq j\leq\Omega.

The synchronized state, where \rho_{j}=1, \phi_{j}=0\quad\forall j, is stable independently of K for 1+c_{1}c_{2}\geq 0, while for 1+c_{1}c_{2}<0 only if the following condition holds:

K>K^{H}_{min}=-\frac{2(1+c_{1}c_{2})}{1+c_{1}^{2}}. | (7) |

Therefore, for each couple (c_{1},c_{2}) we can find a minimum coupling value K^{H}_{min} such that the synchronized state is stable. The resulting stability map is shown in Fig. 4. Notice that the condition 1+c_{1}c_{2}<0 is sufficient for the onset of the instability, when the CGLE is defined on a continuous spatial support [15]. In fact, this is known as the condition of the Benjamin-Feir instability for the CGLE [13, 15].

Stability analysis is more complicated for the splay state. Making use of the recurrence relations (3) we can first compute \rho_{j} and \theta_{j} to evaluate the Jacobian blocks {\bf A}_{j} (see Eqs. (6b)). Then we can assess the stability of the splay state in the plane (c_{1},c_{2}) by computing the Jacobian matrix eigenvalues. An example of the outcome of this procedure is shown in Fig. 5, where the parameters have been set to the values c_{1}=-5, c_{2}=4 and K=4. Here all eigenvalues for j>1 have a negative real part, so that the splay state is linearly stable. Notice the fast convergence of the eigenvalues to their asymptotic state values.

The analysis of the synchronized state and of the splay one of the directed chain of coupled CGL oscilators is summarized in Fig. 6 for the case K=4. The different regions of this diagram are described in the caption; the red cross locates the point in the diagram which defines our working condition as selected in the forthcoming sections when investigating the stochastic version of the directed chain of coupled CGL oscillators. More details on linear stability analysis are given in A.

## 3 Effects of stochasticity

### 3.1 Linear amplification mechanism

The stochastic version of the deterministic model (1) reads

\frac{\mathrm{d}W_{j}}{\mathrm{d}t}\!\!=\!\!W_{j}-(1+ic_{2})|W_{j}|^{2}W_{j}+(% 1+ic_{1})K(W_{j-1}\!-\!W_{j})+\sigma\eta_{j}(t) | (8) |

where \sigma is the noise amplitude, \eta_{j}=\Re(\eta_{j})+i\Im(\eta_{j}) is a complex additive noise with zero mean and correlators \braket{\Re(\eta_{j})(t)\Re(\eta_{l})(t^{\prime})}=\braket{\Im(\eta_{j})(t)\Im% (\eta_{l})(t^{\prime})}=\delta_{jl}\delta(t\!-\!t^{\prime}). In what follows the numerical investigations of the stochastic dynamics (8) has been performed for the parameter values (c_{1},c_{2},K)=(-5,4,4) (see the red cross in Fig. (6) ), where both the synchronized and the splay state of the deterministic dynamics are linearly stable. We want to investigate the effects of a small additive noise on the deterministic evolution (1) [21, 22, 23]. In practice, we have always taken \sigma=10^{-5}, a value which is five orders of magnitude smaller than the oscillations amplitude of the synchronized state. As shown in B, Equation (8) can be rewritten for the polar components of the complex variable W_{j}, while the corresponding noise components remain delta-correlated and – at least near the limit cycle solutions – additive. In practice, we have studied the effects of the noise–induced fluctuations around these states. We know from the previous section that both deterministic states are indeed stable limit cycles with a complex eigenvalues Jacobian. This guarantees the presence of stochastic oscillations, also called quasi-cycles [24, 25], on the top of the deterministic stable states. Then, we can proceed to the Fourier analysis of our system linearized around each limit cycle. We denote by {\bf\delta\tilde{v}} and {\bf\tilde{\xi}} the Fourier transforms of the perturbations vector {\bf\delta v} and of the polar white noise \xi\equiv(\xi_{\rho},\xi_{\theta}), respectively. We can readily obtain {\bf\delta\tilde{v}}_{j}\!\!=\!\!\sum^{2\Omega}_{l=1}\Phi_{jl}^{-1}(\omega){% \bf\tilde{\xi}}_{l}, where \Phi_{jl}\!=\!-J_{jl}-i\omega\delta_{jl}. To pursue the analysis of the oscillations we compute the power spectrum density matrix of the fluctuations in the vicinty of the attractor [26]

\braket{{\bf\delta\tilde{v}}_{l}(\omega){\bf\delta\tilde{v}}_{j}(\omega)}\!\!=% \!\!P_{lj}(\omega)\!\!=\!\!\sum^{2\Omega}_{k=1}\Phi_{lk}^{-1}(\omega)\!\!\left% ({\Phi^{\dagger}}_{kj}\right)^{\!-1}\!\!(\omega). | (9) |

Its diagonal entries are the power spectrum of transversal (j odd) and longitudinal (j even) oscillations around both solutions. We first focus on the transversal, radial, fluctuations around the synchronized state. In Fig. 7(a) we depict the power spectrum of several nodes. The solid line stands for the analytical power spectrum computed from equation (9) while symbols correspond to direct numerical simulations of equation (8), using the Euler-Maruyama algorithm (\mathrm{d}t=0.001). The power spectrum of the first node, peaked at zero frequency (circle, black line) is the one of white noise. As we proceed along the chain, the peak of the power spectrum progrssively shifts towards higher frequencies. The profiles around the peak become narrower (thus singling out a well defined oscillation frequency), while fluctuations are amplified along the chain. This amplification can be well appreciated by direct inspection of Fig. 7(b). Such amplification and modulation proceeds along the chain as long as the linear approximations holds. Out of this approximation, nonlinear effects should take over and stop the amplification process. Note that an analogous phenomenon was already discussed in [19] for noisy fluctuations around a single fixed point. Since the structure of the Jacobian remains essentially the same for the splay state, here we face a qualitatively identical situation. A similar amplification mechanism takes place for longitudinal fluctuations around both stable states, as exemplified in the inset of Fig. 7(a). However, longitudinal oscillations are typically characterized by a broader spectrum, possibly due to the softer nature of the phase direction with respect to the radial one for Ginzburg-Landau potentials. To summarize our findings, noisy fluctuations around both attractors are amplified and modulated as one proceeds along the chain to yield sharper and stronger oscillations. While nonlinear effects would eventually arrest this amplification process, the linear mechanism is typically enough to overcome the attractor linear stability itself. These features are mainly due to the unidirectional structure of the Jacobian, which is highly non-normal. It is well known that non-normality amplifies transient dynamics [27, 28, 29] and may lead to convective instability [20]. Here the presence of noise makes this amplification perpetual [18].

### 3.2 Pattern formation

Why is this so important? Let’s immagine the following scenario where both solutions exist and are stable. We then seed the following initial conditions \rho_{j}(t=0)=1, \phi_{j}(t=0)=0 all over the chain. What we expect from a naive linear stability analysis is that, for small noise amplitudes, the system will remain in the vicinity of the synchronized state, with fluctuations of the order the noise amplitude \sigma. On the contrary, our analysis reveals that the amplification mechanism here discussed will drive the system to progressively explore larger portions of the available phase space, until it eventually reaches the splay state. This is illustrated in Fig. 8, where we show the radial time series of successive nodes. The time series of the first nodes are plotted in red: they remain settled on the synchronized state, the amplification on these first nodes not being strong enough to escape from its basin of attraction. After the 10^{th} node (blue line) fluctuations are now strong enough to escape, and reach the second attractor, settling on the splay state radius \rho_{\bar{j}+1}. Nodes to the right converge to successive radii \rho_{j} with j>\bar{j}. The attractor values \rho_{j}, each represented by a dashed line, is found thanks to the recurence relations (3). They are in good agreement with the time series simulations performed by an Euler-Maruyama algorithm. Obviously, this could not be expected from a traditional linear stability analysis.

Actually, by careful observation of the time series, one can realize that these shifts from the synchronized to the splay state take place as a sort of zipping mechanism from right to left. The rightmost nodes display larger oscillations, and are the first ones to escape the synchronized state (violet lines in Fig. 8). Progressively, the escape point from the synchronized state moves to the left. In this process, individual time series typically jump synchronously from one value \rho_{\bar{j}+k} to the next \rho_{\bar{j}+k+1} (see the green and blue timeseries in Fig. 8).
This zipping process continues towards the left, until the oscillations on node \bar{j} (the 9^{th} node in our case) become strong enough to escape the synchronized part of the attractor, where \rho=1.

A direct consequence of this mechanism is the formation of spatiotemporal patterns [22, 30, 31] as shown in Fig. 9. Our system is initially prepared on the synchronized state and exposed to a noise of amplitude \sigma=10^{-5}. After some time we see that the rightmost nodes easily reach the second attractor. However, as we already discussed, the same amplification and modulation mechanism holds on the splay state. The fluctuations therefore keep on being amplified along the chain allowing the rightmost nodes of our system to travel erratically in phase space. This is exemplified by the blurred part of Fig. 9. Here the mechanism of desynchronization is quite obvious, being the combination of two ingredients: noise and non-normality. While noise is needed to inject some dynamics in the otherwise stable limit cycle, the non-normality is essential to amplify these fluctuations. This is what makes the system deviate from the synchronized to the splay state and then enter an erratic dynamics.

## 4 Conclusion

Noise is often unavoidable and, as such, it should be accommodated for in realistic models of complex natural phenomena. A particularly interesting setting is faced when the stochastic perturbation, being it of endogenous or exogenous origin, resonates with the degree of inherent non-normality. This situation, as displayed by the examined system, yields a self-consistent amplification of the noise component at short times. The resulting growth of the perturbation can in fact drive a symmetry breaking instability, for a choice of the parameters that would instead result in a stable deterministic evolution. In order to dig into this question, we have here examined a directed chain of diffusively coupled, Ginzburg-Landau oscillators. Oscillators are shaked by an external fluctuating drive, of arbitrarily small strength. The system is initiated in a region of parameters where the synchronous solution proves stable, under the deterministic scenario. Working in this setting, we provided analytical and numerical evidence for a noise induced instability which follows the self-consistent amplification of the imposed disturbance across the chain. The limit cycles get modulated along the transversal direction: almost regular, radial oscillations are displayed, which gain in potency node after node. When the transversal modulation gets large enough, oscillators escape the basin of attraction of the synchronized solution, visiting a non trivial attractor, that we have analytically characterized. The interaction between the two attractors yield complex emerging patterns reminiscent of the deterministic Benjamin-Feir instability. The combination of noise and asymmetric couplings can radically alter the limit cycle dynamics: bistability and associated patterns rise, as the noisy signal is dynamically processed, along the unidirectionally coupled chain. It is worth mentioning that an analogous behavior, due to the forward amplification mechanism, is also expected when the (arbitrarely small) noise is only injected in the leftmost node and not on all degrees of freedom as in our current setup. Traditional (deterministic) linear stability analysis is unable to grasp the essence of the phenomenon, an observation which we find particularly relevant given the recent reports on the ubiquity of non-normality in real systems, from communication networks to foodwebs [32]. More refined approaches, such as convective Lyapunov exponents [33, 34, 35, 36] should however be able to predict a convective instability at the purely deterministic level. Resilience to synchronization might prove a valuable asset, exploited to oppose the onset of pathological states, as e.g. epileptic seizures in brain dynamics. Future investigations are planned to shed light onto these families of noise–instigated instabilities, assisted by the non-normal topology of the underlying support, beyond the simplistic case study here considered.

## Acknowledgments

The authors acknowledge financial support from H2020-MSCA-ITN-2015 project COSMOS 642563. We thank Arkady Pikovsky for useful comments.

## Appendix A Linear stability analysis

We now give more details on the linear stability analysis. Consider small perturbations \delta\rho_{j}(t)\ll 1 and \delta\theta_{j}(t)\ll 1 of the limit cycle solutions, W_{j}=(\rho_{j}+\delta\rho_{j})e^{i(\theta_{j}+\delta\theta_{j})}. Linearizing we obtain to first order in the perturbations

\displaystyle\delta\dot{\rho}_{1} | \displaystyle= | \displaystyle-2\delta\rho_{1} | (10) | ||

\displaystyle\delta\dot{\theta}_{1} | \displaystyle= | \displaystyle-2c_{2}\delta\rho_{1} | (11) |

and for j>1

\displaystyle\delta\dot{\rho}_{j} | \displaystyle= | \displaystyle\delta\rho_{j}\left[1-3\rho_{j}^{2}-K\right]+\delta\rho_{j-1}Kf(% \phi_{j})+\left(\delta\theta_{j}-\delta\theta_{j-1}\right)K\rho_{j-1}g(\phi_{j}) | (12) | ||

\displaystyle\delta\dot{\theta}_{j} | \displaystyle= | \displaystyle\delta\rho_{j}\left[-2c_{2}\rho_{j}+K\frac{\rho_{j-1}}{\rho_{j}}g% (\phi_{j})\right]+\delta\rho_{j-1}\frac{K}{\rho_{j}}g(\phi_{j}) | (13) | ||

\displaystyle-\left(\delta\theta_{j}-\delta\theta_{j-1}\right)K\frac{\rho_{j-1% }}{\rho_{j}}f(\phi_{j}) |

where the \rho_{j} and \phi_{j} need to be evaluated on either the
synchronized or the splay state attractor. Obviously zeroth order
terms stemming from the linearization procedure vanish by construction
when evaluated on these two attractors.

Rewriting the linearized equations in a matrix form highlights their simple block structure, due to the unidirectional input from one node to the next. We introduce the 2\Omega dimensional perturbation vector \delta{\bf v}\equiv(\delta\rho_{1},\delta\theta_{1},\delta\rho_{2},\delta% \theta_{2},\ldots,\delta\rho_{\Omega},\delta\theta_{\Omega})^{T} and write

\delta\dot{\bf v}={\bf J}\,\delta{\bf v} | (14) |

where the Jacobian {\bf J} is a 2\Omega\times 2\Omega lower tridiagonal block matrix, composed of 2\times 2 blocks that describe the in-node linearized dynamics ({\bf A} matrices in the following) or the (linearized) interaction with the previous node ({\bf B} matrices).

For instance, in the case of the synchronized state one has

{\bf J}_{H}=\left(\begin{array}[]{ccccc}{\bf A}_{1}&{\bf 0}&{\bf 0}&{\bf 0}&% \ldots\\ {\bf B}_{H}&{\bf A}_{H}&{\bf 0}&{\bf 0}&\ldots\\ {\bf 0}&{\bf B}_{H}&{\bf A}_{H}&{\bf 0}&\ldots\\ {\bf 0}&{\bf 0}&{\bf B}_{H}&{\bf A}_{H}&\ldots\\ \vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right) | (15) |

where

{\bf A}_{1}=\left(\begin{array}[]{lr}-2&0\\ &\\ -2c_{2}&0\end{array}\right) | (16) |

describes the stability of the first uncoupled Landau-Stuart node, while

{\bf A}_{H}=\left(\begin{array}[]{lr}-(2+K)&\;\;Kc_{1}\\ &\\ -(2c_{2}+Kc_{1})&\;\;-K\end{array}\right)\;,\;{\bf B}_{H}=\left(\begin{array}[% ]{lr}K&-Kc_{1}\\ &\\ Kc_{1}&K\end{array}\right) | (17) |

originate from the other nodes (j>1).

Using simple block matrices results one can show that

\det\left({\bf J}_{H}-\lambda{\bf I}_{2\Omega}\right)=\det\left({\bf A}_{1}-% \lambda{\bf I}_{2}\right)\,\left[\det\left({\bf A}_{H}-\lambda{\bf I}_{2}% \right)\right]^{\Omega-1} | (18) |

(where {\bf I}_{h} is the h\times h identity matrix)
so that the eigenvalues of {\bf J}_{H} are given by the ones of {\bf A}_{1} and the ones of {\bf A}_{H} (with multiplicity \Omega-1).

We easily verifiy that {\bf A}_{1} has eigenvalues \lambda_{\rho}=-2 and \lambda_{\theta}=0 and consequently is stable. We are therefore interested in the eigenvalues \lambda_{H} of {\bf A}_{H} that give

\lambda_{H}^{\pm}=-1-K\pm\sqrt{1-2Kc_{1}c_{2}-K^{2}c_{1}^{2}} | (19) |

The real part of the largest eigenvalue \lambda_{H}^{+} has two zeros for K=0 and K=K^{H}_{min} with

K^{H}_{min}=-\frac{2(1+c_{2}c_{1})}{1+c_{1}^{2}} | (20) |

We can determine the stability condition \mathcal{R}\left[\lambda_{H}^{+}\right]<0 by the sign of K^{H}_{min} and of the small K expansion of equation (19),

\lambda_{H}^{+}\approx-K(1+c_{1}c_{2})<0 | (21) |

, which gives the sign of the K derivative of \mathcal{R}\left[\lambda_{H}^{+}\right] near K=0. Note that they are both controlled by the sign of 1+c_{1}c_{2}, so that one immediatly obtains the homogeneous state stability condition given in the main text.

The splay states give rise to a slightly more complicated Jacobian matrices {\bf J}_{S}^{(\bar{j})}. The first \bar{j} rows are identical to the ones of {\bf J}_{H}, while the following ones are obtained evaluating the linearized equation along the splay state part of the attractor. For instance, for \bar{j}=2 we have

{\bf J}_{S}^{(2)}=\left(\begin{array}[]{ccccc}{\bf A}_{1}&{\bf 0}&{\bf 0}&{\bf 0% }&\ldots\\ {\bf B}_{H}&{\bf A}_{H}&{\bf 0}&{\bf 0}&\ldots\\ {\bf 0}&{\bf B}_{3}&{\bf A}_{3}&{\bf 0}&\ldots\\ {\bf 0}&{\bf 0}&{\bf B}_{4}&{\bf A}_{4}&\ldots\\ \vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right) | (22) |

with

{\bf A}_{j}=\left(\begin{array}[]{lr}(1-3\rho_{j}^{2}-K)&\;\;K\rho_{j-1}g(\phi% _{j})\\ &\\ -\left(2c_{2}\rho_{j}+K\frac{\rho_{j-1}}{\rho_{j}}g(\phi_{j})\right)&\;\;-K% \frac{\rho_{j-1}}{\rho_{j}}f(\phi_{j})\end{array}\right) | (23) |

and

{\bf B}_{j}=\left(\begin{array}[]{lr}Kf(\phi_{j})&\;\;\;-K\rho_{j}g(\phi_{j})% \\ &\\ \frac{K}{\rho_{j}}g(\phi_{j})&\;\;\;K\frac{\rho_{j-1}}{\rho_{j}}f(\phi_{j})% \end{array}\right) | (24) |

where functions g and f are defined in the main text in equations (4). Obviusly, for j\gg\bar{j} we have

{\bf A}_{j}\approx{\bf A}_{\infty}=\left(\begin{array}[]{lr}(1-3\rho_{\infty}^% {2}-K)&\;\;K\rho_{\infty}g(\phi_{\infty})\\ &\\ -\left(2c_{2}\rho_{\infty}+Kg(\phi_{\infty})\right)&\;\;-Kf(\phi_{\infty})\end% {array}\right) | (25) |

and

{\bf B}_{j}\approx{\bf B}_{\infty}=\left(\begin{array}[]{lr}Kf(\phi_{\infty})&% \;\;\;-K\rho_{\infty}g(\phi_{\infty})\\ &\\ \frac{K}{\rho_{\infty}}g(\phi_{\infty})&\;\;\;Kf(\phi_{\infty})\end{array}\right) | (26) |

Once again we have

\displaystyle\det\left({\bf J}_{S}^{(\bar{j})}-\lambda{\bf I}_{2\Omega}\right) | \displaystyle\!\!= | \displaystyle\det\left({\bf A}_{1}-\lambda{\bf I}_{2}\right)\,\left[\det\left(% {\bf A}_{H}-\lambda{\bf I}_{2}\right)\right]^{\bar{j}-1} | (27) | ||

\displaystyle\prod_{j=\bar{j}+1}^{N}\left[\det\left({\bf A}_{j}-\lambda{\bf I}% _{2}\right)\right] |

so that to estimate the eigenvalues of {\bf J}_{S}^{(\bar{j})} we also need to compute the eigenvalues of the matrices {\bf A}_{j}, evaluated on the splay attractor values \rho_{j} and \phi_{j} obtained from the recurrence equations (3).

## Appendix B Nature of the noise

In this appendix we shall demonstrate that the additive stochastic corrections we introduced in our system (see equation (8)) remains of the same kind in polar form. We first write the ordinary differential equations for the real and imaginary part of W_{j}=X_{j}+iY_{j}. After few lines of algebra we end up with

\displaystyle\frac{\mathrm{d}X_{j}}{\mathrm{d}t} | \displaystyle\!\!= | \displaystyle X_{j}+\left(X_{j}^{2}+Y_{j}^{2}\right)\left(-X_{j}+c_{2}Y_{j}\right) | (28) | ||

\displaystyle+K\left(X_{j-1}-X_{j}+c_{1}\left(Y_{j}-Y_{j-1}\right)\right) | |||||

\displaystyle+\sigma\eta_{j}^{X} | |||||

\displaystyle\frac{\mathrm{d}Y_{j}}{\mathrm{d}t} | \displaystyle\!\!= | \displaystyle X_{j}+\left(X_{j}^{2}+Y_{j}^{2}\right)\left(-c_{2}X_{j}-Y_{j}\right) | |||

\displaystyle+K\left(Y_{j-1}-Y_{j}+c_{1}\left(X_{j-1}-X_{j}\right)\right) | |||||

\displaystyle+\sigma\eta_{j}^{Y} | |||||

Writing in polar form W_{j}\!\!=\!\!\rho_{j}\exp(i\theta_{j}) implies that \rho_{j}\!\!=\!\!\sqrt{X_{j}^{2}+Y_{j}^{2}} and \theta_{j}\!\!=\!\!Atan\left(\frac{Y_{j}}{X_{j}}\right). In terms of O.D.E.s it means that

\rho_{j}\frac{\mathrm{d}\rho_{j}}{\mathrm{d}t}=X_{j}\frac{\mathrm{d}X_{j}}{% \mathrm{d}t}+Y_{j}\frac{\mathrm{d}Y_{j}}{\mathrm{d}t} | (30a) | ||

\rho_{j}^{2}\frac{\mathrm{d}\theta_{j}}{\mathrm{d}t}=X_{j}\frac{\mathrm{d}Y_{j% }}{\mathrm{d}t}-Y_{j}\frac{\mathrm{d}X_{j}}{\mathrm{d}t} | (30b) |

We now want to obtain the Langevin equation for \rho_{j} and \theta_{j}, this leads to

\displaystyle\rho_{j}\frac{\mathrm{d}\rho_{j}}{\mathrm{d}t} | \displaystyle\!\!= | \displaystyle\rho_{j}^{2}-\rho_{j}^{4}+K\left(-\rho_{j}^{2}+\rho_{j}\rho_{j-1}% f(\phi_{j})\right)+\sigma\left(X_{j}\eta_{j}^{X}+Y_{j}\eta_{j}^{Y}\right) | (31) | ||

\displaystyle\rho_{j}^{2}\frac{\mathrm{d}\theta_{j}}{\mathrm{d}t} | \displaystyle\!\!= | \displaystyle-c_{2}\rho_{j}^{2}+K\left(-c_{1}\rho_{j}^{2}+\rho_{j}\rho_{j-1}g(% \phi_{j})\right)+\sigma\left(X_{j}\eta_{j}^{Y}-Y_{j}\eta_{j}^{X}\right) | |||

where the auxiliary functions f and g have been introduced in equations (4). The sum of two Gaussian variable is itself a Gaussian variable, whose average value is the sum of the two previous average value while its variance is the quadratic sum of the variances. Therefore we can introduce two new Gaussian delta correlated and zero mean white noise variables \xi_{j}^{\rho} and \xi_{j}^{\theta} such that their standard deviations are

\Sigma_{\rho,\theta}=\sqrt{X_{j}^{2}+Y_{j}^{2}}=\rho_{j} | (33) |

This leads to the final Langevin equations in polar form

\displaystyle\frac{\mathrm{d}\rho_{j}}{\mathrm{d}t} | \displaystyle\!\!= | \displaystyle\rho_{j}-\rho_{j}^{3}+K\left(-\rho_{j}+\rho_{j-1}f(\phi_{j})% \right)+\sigma\xi^{\rho}_{j} | (34) | ||

\displaystyle\frac{\mathrm{d}\theta_{j}}{\mathrm{d}t} | \displaystyle\!\!= | \displaystyle-c_{2}+K\left(-c_{1}+\frac{\rho_{j-1}}{\rho_{j}}g(\phi_{j})\right% )+\frac{\sigma}{\rho_{j}}\xi^{\theta}_{j} | |||

which display a multiplicative but delta correlated zero-average noisy term. In our power spectrum analysis, conducted expanding near the limit cycle solutions, this multiplicative component can be safely approximated by its limit cycle value, making the dominant noise component additive.

## References

- [1] Kuramoto Y 1984 Chemical Oscillations, Waves, and Turbulence Springer-Verlag New York
- [2] Nicolis G and Prigogine I, Self-Organization in Non-Equilibrium Systems 1977 Wiley
- [3] Goldbeter Biochemical Oscillations and Cellular Rhythms 1997 (Cambridge: University Press Cambridge)
- [4] Pikovsky A, Rosenblum M, and Kurths J 2001 Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge: Cambridge University Press)
- [5] Latora V, Nicosia V, and Russo G 2017 Complex Networks: Principles, Methods and Applications (Cambridge: Cambridge University Press)
- [6] Barrat A, Barthelemy M, Vespignani A 2008 Dynamical processes on complex networks (Cambridge: Cambridge University Press)
- [7] Cross M C and Hohenberg P C 1993 Pattern formation outside of equilibrium Rev. Mod. Phys. 65 851
- [8] Nakao H 2014 Eur. Phys. J. Special Topics 223 2411
- [9] Pecora L M and Carroll T L 1998 Master Stability Functions for Synchronized Coupled Systems Phys. Rev. Lett. 80 2109
- [10] Koseska A, Volkov E and Kurths J 2013 Oscillation quenching mechanisms: Amplitude vs Oscillation Death Phys. Rep. 531 173
- [11] Challenger J D, Burioni R, and Fanelli D Turing-like instabilities from a limit cycle Phys. Rev. E 92 022818
- [12] Aranson I S and Kramer L 2002 The world of the complex Ginzburg-Landau equation Rev. Mod. Phys. 74 99
- [13] Benjamin T B and Feir J E 1967 The disintegration of wave trains on deep water J. Fluid. Mech. 27 417
- [14] Stuart J T and DiPrima R C 1978 The Eckhaus and Benjamin-Feir resonance mechanisms Proc. R. Soc. Lond. A. 362 27
- [15] Di Patti F, Fanelli D, Carletti T 2016 Drift-induced Benjamin-Feir instabilities Europhys. Lett. 114 6
- [16] Di Patti F, Fanelli D, Miele F and Carletti T 2017 Benjamin-Feir instabilities on directed networks Chaos, Solitons and Fractals 96 8-16
- [17] Trefethen L N and Embree M 2005 Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators (Princeton: Princeton University Press)
- [18] Nicoletti S, Zagli N, Fanelli D, Livi R, Carletti T and Innocenti G 2018 Non-normal amplification of stochastic quasicycles Phys. Rev. E. in press .
- [19] Fanelli D, Ginelli F, Livi R, Zagli N and Zankoc C 2018 Noise driven neuromorphic tuned amplifier Phys. Rev. E 96 062313
- [20] Deissler R J 1989 External noise and the origin and dynamics of structure in convectively unstable systems Journal of Statistical Physics 54 1459-1488
- [21] Asllani M, Biancalani T, Fanelli D and McKane A J 2013 The linear noise approximation for reaction-diffusion systems on networks The European Physical Journal B
- [22] McKane A J, Biancalani T and Rogers T 2014 Stochastic pattern formation and spontaneous polarisation: the linear noise approximation and beyond Bull Math Biol. 76(4):895-921.
- [23] Boland R P, Galla T and McKane A J 2008 How limit cycles and quasi-cycles are related in systems with intrinsic noise Journal of Statistical Mechanics: Theory and Experiment 2008
- [24] Lugo C A and McKane A J 2008 Quasicycles in a spatial predator-prey model Phys. Rev. E 78 051911 (2008)
- [25] McKane A J and Newman T J 2005 Predator-Prey Cycles from Resonant Amplification of Demographic Stochasticity Phys. Rev. Lett. 94 218102
- [26] Challenger J D and McKane A J 2013Synchronization of stochastic oscillators in biochemical systems Phys. Rev. E 88 012107
- [27] Neubert M G and Caswell H 1997 Alternatives to resilience for measuring the responses of ecological systems to perturbations Ecology 78 653-665.
- [28] Trefethen L N, Trefethen A E, Reddy S C and Driscoll T A 1993 Science 261 5121 578-584
- [29] Lewis M A and Kareiva P 1993 Allee Dynamics and the Spread of Invading Organisms Theoretical Population Biology 43 2 141-158
- [30] Asslani M, Di Patti F and Fanelli D 2012 Stochastic Turing patterns on a network Phys. Rev. E 86 046105
- [31] Biancalani T, Jafarpour F and Goldenfeld N 2017 Giant Amplification of Noise in Fluctuation-Induced Pattern Formation Phys. Rev. Lett. 118 018101
- [32] Asllani M and Carletti T 2018 Universality of non-normality in real complex networks, ArXiv, 1803.11542
- [33] Deissler R J and Kaneko K 1987 Phys Lett. 119 397 (1987)
- [34] Lepri S, Politi A and Torcini A 1996 J. Stat. Phys. 82 1429
- [35] Lepri S, Politi A and Torcini A 1997 J. Stat. Phys. 88 31
- [36] Kenfack Jiotsa A, Politi A and Torcini A 2013 Convective Lyapunov spectra Journal of Physics A: Mathematical and Theoretical 46 25 254013