We are interested in the situation where, due to an external perturbation or change in the environment, the intrinsic timescales of the individual oscillators separate into a fast and a slow subsystem. We consider a set of N oscillators each with a compact phase degree of freedom θ _i ∈ (−π, π] whose time-evolution is governed by the set of coupled differential equations (Kuramoto, 1975), d i θ ̇ i = ω i − ∑ j = 1 N b i j sin θ i − θ j + η i , (1) for i = 1, … N. The natural frequency of the ith oscillator is denoted ω _i , the structure of the coupling network is given by elements b _ij of the adjacency matrix (Newman, 2018). Ambient noise is modelled at the ith oscillator by η _i and is taken as a time-correlated noise, uncorrelated in space, i.e., ⟨ η i ( t ) η j ( t ′ ) ⟩ = η 0 , i 2 δ i j ⁡ exp [ − | t − t ′ | / τ i ] , where τ _i is the typical correlation time of the noise at the ith oscillator. The non-negative parameters d _i ’s define the individual timescale of each oscillator. Removing the noise term, Eq. 1 may have multiple stable fixed points of the dynamics which essentially depend on the coupling topology and strength as well as the distribution of natural frequencies. Below, we assume that such a stable fixed point { θ i ∗ } exists and that the noise term is small enough such that the dynamics remains inside the initial basin of attraction. In the present scenario, we assume that we have two sets of oscillators that we denote F and S , respectively with N F and N S oscillators, such that d i = d ̄ i ∈ F d ̲ i ∈ S (2) with d ̄ ≪ d ̲ . The latter means that oscillators belonging to S have a much slower intrinsic timescale than those belonging to F . In the following, we focus on the dynamics of the oscillators in the slow component. Within the assumption of timescale separation, one can rewrite Eq. 1 as, θ ̇ i = ω i − ∑ j = 1 N b i j sin θ i − θ j + η i , i ∈ S ϵ θ ̇ i = ω i − ∑ j = 1 N b i j sin θ i − θ j + η i , i ∈ F , (3)

where we defined d ̲ / d ̄ = ϵ − 1 and, without loss of generality, set d ̲ to unity. In the limit ϵ → 0, the oscillators within F instantaneously adapt their phases. In the next section, we consider the dynamical system Eq. 3 in the vicinity of a stable fixed point and perform a singular perturbation analysis using Mori-Zwanzig formalism.

Even though we consider Kuramoto oscillators, the following approach applies in general to coupled dynamical systems that have a stable fixed point around which they evolve and where linearization is valid. To analyze the resilience of the slow component, we consider the dynamics of the system close to a fixed point { θ i ∗ } . In particular, we are interested in the time-evolution of the phase deviations x i ( t ) = θ i ( t ) − θ i ∗ for i ∈ S and y i ( t ) = θ i ( t ) − θ i ∗ for i ∈ F whose dynamics at the first order reads, x ̇ ϵ y ̇ = J S S J S F J F S J F F x y + η S η F (4) where we defined the matrix J i j = b i j ⁡ cos θ i ∗ − θ j ∗ i ≠ j − ∑ k = 1 N b i k ⁡ cos θ i ∗ − θ k ∗ i = j , (5) which is the Jacobian of the system and is a Laplacian matrix when phase differences are between − π 2 and π 2 . Using Mori-Zwanzig formalism (Mori, 1965; Zwanzig, 1973) with x and y being respectively the resolved and unresolved variables (see Supplementary Material Appendix A), one can express the first row of Eq. 4 as, x ̇ i = ∑ j = 1 N S J S S i j x j + η S i + ∑ α = 1 N F ∫ 0 t ϵ − 1 e ν α ϵ − 1 t − t ′ × ∑ k , m = 1 N F ∑ l = 1 N S J F S k l x l t ′ w α , k J S F i m w α , m d t ′ + ∑ α = 1 N F ∫ 0 t ϵ − 1 e ν α ϵ − 1 t − t ′ ∑ k , m = 1 N F η F k t ′ w α , k J S F i m w α , m d t ′ , (6) where we denoted w _α the eigenvectors of J F F with corresponding eigenvalues ν _α < 0. In Eq. 6, the first two terms on the right-hand side are Markovian, i.e., they do not depend on the history of the trajectories, while the other ones have memory and thus do depend on the history. We are interested in the time-evolution of the slow component x when there is a timescale separation, i.e., ϵ → 0. Taking the latter limit in Eq. 6, and using lim φ → ∞ ∫ 0 t φ e φ t ′ f ( t ′ ) d t ′ = f ( 0 ) (see Supplementary Material Appendix B), and the identity J F F − 1 = ∑ α ν α − 1 w α w α ⊤ , yields in a matrix form, x ̇ = J S S x − J S F J F F − 1 J F S x + η S − J S F J F F − 1 η F = J red x + ξ , (7) where in the second line we defined the reduced Jacobian J red = J S S − J S F J F F − 1 J F S , and denoted the noise term as ξ = η S − J S F J F F − 1 η F . It is interesting to note that the reduced dynamics given by Eq. 7 can be obtained by a Kron reduction (Kron, 1939; Dörfler and Bullo, 2012) of the fast component of the system (see (Tyloo et al., 2023) for an example). The dynamics of the slow component is then governed by Eq. 7 where the effective noise at the ith oscillator is a combination of the noise at the ith oscillator with a superposition of the noise inputs at oscillators belonging to the fast component. For undirected coupling as we consider in the following, one has that J F S = J S F ⊤ . The linear system Eq. 7 can be solved by expanding over the eigenmodes of J _red denoted u _α , with corresponding eigenvalues λ _α , α = 1 , … N S , i.e., x(t) = ∑ _α c _α (t)u _α . As J _red is also the negative of a Laplacian matrix, one has that 0 = λ 1 > … ≥ λ N S with u 1 , i = 1 / N S . The expansion coefficients satisfy the uncoupled differential equations, c ̇ α = λ α c α + ξ ⋅ u α , α = 1 , … , N S . (8) Assuming a vanishing initial condition, the general solution to Eq. 7 is given by, x i t = ∑ α = 2 N S ∫ 0 t ⁡ exp λ α t − t ′ u α ⋅ ξ t ′ d t ′ u α , i + t N S ∑ j = 1 N S ξ j , (9) for i = 1 , … N S . One remarks that, if ∑ j = 1 N S ξ j ≠ 0 , then all the oscillators in the slow component drift simultaneously together in time, i.e., along u ₁. Such a homogeneous overall shift in the phases of all the oscillators in the slow component does not change the system. Indeed, this is due to the global rotational symmetry of the original system Eq. 3, and the timescale separation. Accordingly, in the following, all the averages over time or noise sequences, denoted ⟨.⟩, correspond to contributions orthogonal to u ₁. The expression Eq. 9 can be used directly to calculate the moments of the phase deviations.

Various characteristics of the response can be used to determine the resilience of the coupled oscillators. When subject to stochastic inputs, a natural choice is to evaluate the magnitude of the deviations from the synchronized fixed point by calculating the variance of the phase deviations. Here, we consider time-correlated noise of the form given below Eq. 1, with τ i = τ S if i ∈ S and τ i = τ F if i ∈ F , where τ S , τ F are the typical correlation times respectively, in the slow and fast component. This noise translates into the reduced noise as ⟨ ξ j ( t 1 ) ξ k ( t 2 ) ⟩ = η 0 2 δ j k ⁡ exp ( − | t 1 − t 2 | / τ S ) + η 0 2 [ J S F J F F − 2 J F S ] j k ⁡ exp ( − | t 1 − t 2 | / τ F ) . The variance of the phase deviations in the slow component is defined as ⟨ ( x i − N S − 1 ∑ j ∈ S x j ) 2 ⟩ and simply denoted as ⟨ x i 2 ⟩ in the rest of this manuscript to shorten the notation. It is calculated from Eq. 9 and reads in the long time limit (see Supplementary Material Appendix C), ⟨ x i 2 ⟩ = η 0 , S 2 ∑ α = 2 N S u α , i 2 λ α λ α − τ S − 1 + η 0 , F 2 ∑ α , β = 2 N S λ α + λ β − 2 τ F − 1 Γ α β u α , i u β , i λ α + λ β τ F − 1 − λ α τ F − 1 − λ β , (10) with the scalar Γ α β = u α ⊤ J S F J F F − 2 J F S u β . In Eq. 10, we set the standard deviation of the ambient noise in the slow and fast components respectively, to η 0 , S and η 0 , F . We also set distinct homogeneous correlation times for the noise in each component as τ i = τ S for i ∈ S and τ i = τ F for i ∈ F . While the contribution to the variance from the additive noise in the slow component is essentially given by the position of the oscillators on the slowest eigenmodes, the effect of the noise coming from the fast component involves combinations of eigenmodes. The precise combination depends on the effective reduced dynamics through Γ _αβ . The shortest timescale in the system is set by the oscillators belonging to F . However, by tuning the correlation time of the noise τ = τ F = τ S , one can investigate the regimes where the λ N S τ ≪ 1 and λ ₂ τ ≫ 1. Indeed, in the limit where the noise correlation time is shorter than the timescales of the slow component and the same in both components, the variance becomes, ⟨ x i 2 ⟩ = η 0 , S 2 τ ∑ α = 2 N S u α , i 2 − λ α + η 0 , F 2 τ ∑ α , β = 2 N S 2 Γ α β u α , i u β , i − λ α − λ β . (11) In the other limit where the noise correlation time is the longest timescale, one has, ⟨ x i 2 ⟩ = η 0 , S 2 ∑ α = 2 N S u α , i 2 λ α 2 + η 0 , F 2 ∑ α , β = 2 N S Γ α β u α , i u β , i λ α λ β . (12)

Comparing the two limiting cases Eqs 11, 12, one remarks that in both variances, a significant contribution might come from the slowest eigenmodes. Note also that Eq. 10 is more generally valid in the case where τ F and τ S are different.

In the Supplementary Material Appendix D, we give the variance of the phases when there is no timescale separation.

To obtain more insights into the contribution from the fast component, let us consider specific situations in the strong coupling limit.

In the strong coupling limit, one has | θ i * − θ j * | ≪ 1 ∀ i , j , so that one can approximate the Jacobian Eq. 5 as, J i j = b i j i ≠ j − ∑ k = 1 N b i k i = j . (13)

Within this coupling limit and some other assumptions that are specified below, one can further consider network structures that give more insights about Eq. 10.

We start by checking the analytical prediction when the oscillators in both components are randomly chosen. In Figure 5, we consider a Watts-Strogatz network (Newman, 2018) with m = 4 initial nearest neighbors and a rewiring probability p _rewiring = 0.1, of size N S + N F = 40 , with respectively N S = 21 and N F = 19 oscillators in each component. The timescale separation is numerically simulated by taking the limit ϵ = d ̄ / d ̲ ≪ 1 . In Figure 5, we show the variance of the oscillators belonging to S for different values of 0 < ϵ ≤ 1. One observes that, for ϵ = 1, i.e., no timescale separation such that all the oscillators belong to S , the variance follows the analytical prediction of Eq. D1 in the Supplementary Material. Note that we are only showing the oscillators that belongs the slow component when ϵ < 1. When ϵ = 0.5, the variance differ from both prediction Eqs. (D1), (10). Eventually, for ϵ = 0.01, the numerical simulations follow the analytical prediction of Eq. 10, which corresponds to a system with a timescale separation. Besides confirming the theoretical prediction, Figure 5 shows that a timescale separation can induce significant changes in the variance of the oscillators in the slow component. Indeed, the variance of oscillator 19 when ϵ = 1 is approximately three times smaller than when ϵ = 0.01.

We then move to the case described Section 2.4 where the oscillators in the fast component are disconnected. In Figure 6, we numerically simulate the dynamics of Eq. 1 on a Watts-Strogatz network with m = 4 initial nearest neighbors and a rewiring probability p _rewiring = 0.3, of size N S + N F = 20 , with respectively N S = 17 and N F = 3 oscillators in each component. For the oscillators in the fast component, we picked the two oscillators with highest and lowest degrees, plus an additional one with an intermediate degree. As expected, we find that the effect of the timescale separation is limited to the oscillators directly connected to those in the fast compoenent. Overall, compared to the case without any timescale separation, the variance is not much affected.

We consider the situation described in Section 2.4.2 where oscillators in the fast component that have some connections to the slow ones, are connected to all of them. The numerical results are shown on Figure 7 where this particular setting has been simulated for a modified Erdős-Rényi network of N S + N F = 13 oscillators with various correlation times of the noise and different levels of heterogeneity in the natural frequencies. Here, only the fast component is subjected to noise, while the slow component is noiseless. One observes that increasing the heterogeneity in the natural frequencies of the oscillators induces larger variances for the phase deviations. Not shown on Figure 7 is the homogeneous case of oscillators with identical natural frequencies, for which the variance vanishes as predicted in Section 2.4.2. While the heterogeneity increases the noise transmission from the fast to the slow component, one observes that the amplitude of the deviations is still rather small in Figure 7. Besides u α ⊤ J S F being small, this is because the oscillators in the fast component with inter-component connections have relatively large degrees, which directly reduces the noise transmission in Γ α β = u α ⊤ J S F J F F − 2 J F S u β .

In the other situation where some oscillators in the slow components are connected to a large fraction of the oscillators in the fast component, we showed in Section 2.4.3 that their variance is more important than oscillators with fewer or no connection to the fast component. This result is particularly interesting and intriguing, as in the regular situation where there is no timescale separation, oscillators that have a larger number of connections to the other elements typically have a smaller variance (Tyloo et al., 2019). Indeed, this is first illustrated in Figure 8 where the variance of each oscillator is given by the color map when there is no timescale separation (left top panels), and when there is a timescale separation (right top panels). In both correlation time limits (with τ S = τ F ), one observes that oscillators belonging to S having smaller variances in the left top and middle panels are among the ones with larger variances in the right top and middle panels. In particular for the long correlation time limit, one also observes that the timescale separation modifies the variance of the oscillators far from the fast component. The theory and the numerical simulations for the two systems with and without timescale separation are confirmed in the bottom panel of Figure 8. We then move to the case where the noise correlation times are different in each component. In Figure 9, we show both limits λ N S τ F ≫ 1 and λ 2 τ S ≪ 1 (top panels), and λ N S τ F ≪ 1 and λ 2 τ S ≫ 1 (middle panels). Similar amplification of the fluctuations as in the previous case are observed. However, comparing the left and right middle panels, one remarks that some oscillators well connected to the fast component keep rather small variances while others having fewer connections become more vulnerable.

Here, we check the theory when the system is defined on a layered network, i.e., the slow and fast components each corresponds to one layer. We consider the specific setting where both layers have the same network connectivity and the inter-layer coupling is made through single connections between corresponding oscillators in each component (such structure are sometimes called multiplex (Newman, 2018)). In Figure 10, the numerical simulations for the variance (dots and crosses) match the theory Eqs. 10 and Eq. D1 in the Supplementary Material (solid and dashed lines) for various correlation times of the noise (homogeneous, i.e., τ F = τ S ). The blue and orange data points correspond to the situation where there is a timescale separation between the two components, while for the red and green ones, there is no timescale separation. Interestingly, one observes that the two different situations produce similar variances for the phases. As predicted in Section 2.4.4, in the limit where τ is the longest timescale in the system, the variances in the two situations are close to each other (green and red data points in Figure 10). A similar effect is observed in the short correlation time limit, where the variance is only marginally impacted by the timescale separation. We do not show but found similar behavior in the case of heterogeneous noise correaltion time in the two components. This indicates that networks having a layered structure are robust to timescale separation when the slow and fast components are defined on each layer.