Lyapunov Analysis of Chimera States in Globally Coupled Stuart-Landau oscillators

TL;DR

This paper uses Lyapunov analysis to explore the collective dynamics of chimera states in globally coupled Stuart-Landau oscillators, revealing differences in their Lyapunov spectra and modes depending on the type of global coupling.

Contribution

It introduces a Lyapunov-based approach to distinguish and analyze chimera states in different globally coupled oscillator systems, highlighting qualitative differences in their collective modes.

Findings

Lyapunov spectra split into distinct groups with system size

Lyapunov dimension scales extensively with system size

Different coupling types exhibit distinct collective Lyapunov modes

Abstract

Oscillatory systems with long-range or global coupling offer promising insight into the interplay between high-dimensional (or microscopic) chaotic motion and collective interaction patterns. Within this paper, we use Lyapunov analysis to investigate whether chimera states in globally coupled Stuart-Landau (SL) oscillators exhibit collective degrees of freedom. We compare two types of chimera states, which emerge in SL ensembles with linear and nonlinear global coupling, respectively, the latter introducing a constraint that conserves the oscillation of the mean. Lyapunov spectra reveal that for both chimera states the Lyapunov exponents split into different groups with different convergence properties in the limit of large system size. Furthermore, in both cases the Lyapunov dimension is found to scale extensively and the localization properties of covariant Lypunov vectors manifest the presence of collective Lyapunov modes. Here, however, we find qualitative differences between the two types of chimera states: Whereas the ones in the system under nonlinear global coupling exhibit only slow collective modes corresponding to Lyapunov exponents equal or close to zero, those which experience the linear mean-field coupling exhibit also faster collective modes associated with Lyapunov exponents with large positive or negative values.

Figures

Equations

$${\partial }_t{W}_k=W_k-\left(1+i{c}_2\right){\left|W_k\right|}^2{W}_k+\kappa \left(\frac{1}{N}\sum _{j=1}^N{W}_j-W_k\right),$$

$${\partial }_t{W}_k$$

$$-\left(1+i{c}_2\right)\left(\frac{1}{N}\sum _{j=1}^N{\left|W_j\right|}^2{W}_j-{\left|W_k\right|}^2{W}_k\right),$$

$$\left\lvert\left\langle W\right\rangle\right\rvert=\bigg{\lvert}\sum_{j}W_{k}(t=0)/N\bigg{\rvert}=\eta=0.67\text{.}$$

$$\dot{\mathbf{x}}=\mathbf{f}\left(\mathbf{x}\right)\text{.}$$

$$a_k=\frac{1}{\sqrt{2}}(W_k+W_k^{\ast })\text{,}$$

$$b_k=\frac{1}{\sqrt{2}}(W_k-W_k^{\ast })\text{,}$$

$$\bm{x}=\left(\begin{array}[]{c}\bm{a}\\ \bm{b}\end{array}\right)=\left(a_{1},\ldots,a_{N},b_{1},\ldots,b_{N}\right)^{T}\text{.}$$

$$\delta \dot{\mathbf{x}}={\frac{\partial \mathbf{f}}{\partial \mathbf{x}}|}_{\hat{\mathbf{x}}\left(t\right)}\delta \mathbf{x}\text{,}$$

$$\delta \mathbf{x}\left(t\right)={\Psi }_{\hat{\mathbf{x}}}\left(t,t_0\right)\delta {\mathbf{x}}_0\text{,}$$

$${\Vert \delta \mathbf{x}\left(t\right)\Vert }^2=\delta {\mathbf{x}}_0^T{\Psi }_{\hat{\mathbf{x}}}{\left(t,t_0\right)}^T{\Psi }_{\hat{\mathbf{x}}}\left(t,t_0\right)\delta {\mathbf{x}}_0\text{.}$$

$$M_{\hat{\mathbf{x}}\left(t_0\right)}^{(\pm )}=\underset{t\to \pm \infty }{\lim }{\left[{\Psi }_{\hat{\mathbf{x}}}{\left(t,t_0\right)}^T{\Psi }_{\hat{\mathbf{x}}}\left(t,t_0\right)\right]}^{\frac{1}{2t}}\text{,}$$

$$\Vert {\Psi }_{\hat{\mathbf{x}}}\left(t,t_0\right)\delta {\mathbf{x}}_0\Vert \sim e^{{\lambda }_{k}t}\Vert \delta {\mathbf{x}}_0\Vert \text{,}$$

$${\Gamma }_{\hat{\mathbf{x}}\left(t_0\right),k}^{(+)}$$

$${\Gamma }_{\hat{\mathbf{x}}\left(t_0\right),k}^{(-)}$$

$${\Omega }_{\hat{\mathbf{x}}\left(t_0\right),k}={\Gamma }_{\hat{\mathbf{x}}\left(t_0\right),k}^{(+)}\cap {\Gamma }_{\hat{\mathbf{x}}\left(t_0\right),k}^{(-)}\text{.}$$

$$T_{\hat{\mathbf{x}}\left(t_0\right)}{\mathbb{R}}^{2N}=\underset{k=1}{\overset{r}{⨁}}{\Omega }_{\hat{\mathbf{x}}\left(t_0\right),k}\text{.}$$

$${\Omega }_{\hat{\mathbf{x}}\left(t\right),k}={\Psi }_{\hat{\mathbf{x}}}\left(t,t_0\right){\Omega }_{\hat{\mathbf{x}}\left(t_0\right),k}\text{.}$$

$$D_{\text{L}}:=L+\frac{{\sum }_{l=1}^L{\lambda }_l}{\left|{\lambda }_{L+1}\right|}\text{.}$$

$$\delta \mathbf{x}={\left(\delta a_1,\dots ,\delta a_N,\delta b_1,\dots ,\delta b_N\right)}^T\in {\mathbb{R}}^{2N}\text{,}$$

$$\text{IPR}=\sum _{k=1}^N{\left(\delta a_k^2+\delta b_k^2\right)}^2\text{.}$$

$$\tilde{l}:=(l-1/2)/(2N)\in \left(0,1\right)\text{,}$$

$$\delta {\mathbf{x}}_{\text{ts}}\propto \mathbf{f}\left(\hat{\mathbf{x}}\left(t\right)\right)\text{,}$$

$$\delta {\mathbf{x}}_{\text{ps}}\propto {\left(-b_1,\dots ,-b_N,a_1,\dots ,a_N\right)}^T\text{,}$$

$$d/dt\left|⟨W⟩\right|=0\text{.}$$

$$d⟨W⟩/dt=-\nu ⟨W⟩\text{,}$$

$$\dot{\bm{x}}=\left(\begin{array}[]{c}\dot{\bm{a}}\\ \dot{\bm{b}}\end{array}\right)=\left(\begin{array}[]{c}\bm{f}^{(a)}\!\left(\bm{a},\bm{b}\right)\\ \bm{f}^{(b)}\!\left(\bm{a},\bm{b}\right)\text{.}\end{array}\right)=\bm{f}\!\left(\bm{x}\right)\text{,}$$

$$\dfrac{\partial\bm{f}}{\partial\bm{x}}=\left(\begin{array}[]{c | c}\dfrac{\partial\bm{f}^{(a)}}{\partial\bm{a}}&\dfrac{\partial\bm{f}^{(a)}}{\partial\bm{b}}\\ \hline\cr\dfrac{\partial\bm{f}^{(b)}}{\partial\bm{a}}&\dfrac{\partial\bm{f}^{(b)}}{\partial\bm{b}}\end{array}\right)\text{.}$$

$$f_k^{(a)}\left(\mathbf{a},\mathbf{b}\right)=g^{(a)}\left(a_k,b_k\right)+\sum _{j=1}^N{c}^{(a)}\left(a_j,b_j\right)\text{,}$$

$$\frac{\partial\bm{f}^{(p)}}{\partial\bm{q}}=\left(\begin{array}[]{c | c}g_{pq}\bm{1}+c_{pq}E&\dots\\ \hline\cr R_{pq}&\dots\end{array}\right)\text{,}$$

$$g_{pq}$$

$$c_{pq}$$

Full text

Lyapunov exponents; covariant Lyapunov vector; Lyapunov dimension; Lyapunov analysis; chimera states; Stuart-Landau oscillator; global coupling

I Introduction

Chaos readily emerges in systems composed of many coupled oscillatory ’units’. These units may represent individual oscillators or infinitesimally small patches of spatially extended oscillatory media, and the coupling is what gives rise to chaotic instabilities. The origin of the instabilities, however, might lie in microscopic, i.e. local, interactions, involving only a small number of oscillators, or may result from complex interaction patterns on macroscopic system scales. The character of the corresponding chaotic dynamics, accordingly, may vary significantly.

One way to shine light on the nature of the chaotic dynamics is to study chaotic states in oscillatory ensembles containing different numbers of oscillators $N$ while being subject to equivalent coupling schemes. Thereby, the chaotic dynamics are usually characterized in terms of Lyapunov exponents (LEs), which measure the average infinitesimal divergence rates of the motion in phase space Oseledets (1968); Pikovsky and Politi (2016). Our understanding of chaos within this context is comparatively advanced in two limiting cases, so-called intensive and solely extensive chaos. Intensive chaos is characterized by an (in general small) number of positive LEs that is independent of the number of oscillatory units. It occurs, e.g., for certain parameters in systems of globally coupled oscillators  Matthews et al. (1991); Hakim and Rappel (1992); Nakagawa and Kuramoto (1993). In the opposite case, when the chaotic motion arises exclusively from microscopic degrees of freedom and the number of positive LEs scales linearly with $N$, then the chaotic dynamics is solely extensive. Such behavior was shown to exist in some generic spatially one-dimensional models Livi et al. (1986). However, it has been recently established that large systems with global couplings schemes or spatially extended systems in two or three spatial dimensions, might exhibit chaotic motion that belongs to neither of these two situations. Instead, their dynamics is characterized by both a few, collective macroscopic modes and a large number of microscopic, chaotic degrees of freedom Matthews et al. (1991); Hakim and Rappel (1992); Nakagawa and Kuramoto (1993). Our knowledge how and when such collective modes develop and which role they play in the overall high-dimensional chaotic dynamics is still very limited.

The investigation of collective dynamics in high-dimensional chaotic systems remains a challenge. For a long time, it was common belief that conventional Lyapunov analysis cannot capture these collectively chaotic dynamics correctly Shibata and Kaneko (1998). Instead, finite-size perturbations needed to be studied in order to identify collective dynamical modes and the associated LEs. Only recently, Takeuchi et al. Takeuchi et al. (2009); Takeuchi and Chaté (2013) were able to demonstrate that standard Lyapunov analysis can in fact provide information on collective dynamics. The key of their analysis was the calculation of the covariant Lyapunov vectors (CLVs) associated to the LEs Oseledets (1968); Ginelli et al. (2007); Wolfe and Samelson (2007). Having access to CLVs, Takeuchi et al. complemented the information provided by the Lyapunov spectrum with an investigation of the localization or delocalization properties of CLVs associated with particular LEs. They illustrated their strategy for $N$ globally coupled Stuart-Landau (SL) oscillators in a high-dimensional chaotic state. Most of the CLVs were found to be well localized for a large range of ensemble sizes $N$, while a small number of modes appeared to become increasingly delocalized with increasing $N$. The respective perturbations, named collective Lyapunov modes, were shown to be related to macroscopically chaotic degrees of freedom. Moreover, in accordance with earlier work on collective chaos by Nakagawa and Kuramoto Nakagawa and Kuramoto (1993, 1994, 1995), they showed that the overall Lyapunov spectrum can be separated into parts with different convergence properties in the limit of large $N$. In particular, they identified an extensively scaling group of $\mathcal{O}\!\left(N\right)$ positive LEs in the middle of the spectrum, resulting in a flattening of the spectrum. Around the most positive and most negative exponents, however, they also observed sub-extensive groups of exponents, the number of which was shown to scale approximately as $\mathcal{O}\!\left(\ln\!\left(N\right)\right)$ Takeuchi and Chaté (2013).

In this paper, we perform a related Lyapunov analysis of chimera states in globally coupled SL oscillators. Being composed of coexisting groups exhibiting synchronous and asynchronous motion, respectively, chimera states are peculiar dynamic states which can be seen as ’a natural link between coherence and incoherence’ Omel’chenko et al. (2008). They might offer insight into natural phenomena such as some neural activity patterns Sakaguchi (2006) or hydrodynamic flows with mixed laminar and turbulent patterns Barkley and Tuckerman (2005). So far, there are only a few studies on Lyapunov analysis of chimera states that mainly focus on systems of coupled phase oscillators. Wolfrum et al. Wolfrum et al. (2011); Wolfrum and Omel’chenko (2011) showed that the Lyapunov spectra for chimera states in coupled phase oscillator networks of finite size exhibit an extensive number of positive LEs, revealing the hyperchaotic nature of chimera states. Yet, all the positive exponents decay with increasing system size and finally yield neutrally stable zero exponents in the limit $N\to\infty$. Nevertheless, they found the Lyapunov dimension $D_{\text{L}}$ to scale extensively with system size. Furthermore, Botha et al. Botha (2016); Botha and Kolahchi (2018) studied the distribution of finite-time LEs in chimera-like dynamics and identified characteristic patterns in their temporal distribution function. While for usual chaotic states without synchronization pattern the finite-time LEs follow a Gaussian distribution, they found that the distribution in the case of chimera states possesses a complex, multi-modal shape, which they proposed to use as an indicator for chimera-like behavior.

In this work, however, we want to focus our attention on properties of the asymptotic Lyapunov spectra, without further reference to the temporal distribution of finite-time exponents. In particular, we examine whether Lyapunov spectra for chimera states of the globally coupled SL ensemble exhibit similar patterns as observed for the above-mentioned chaotic states, which exist in the same type of oscillatory network for different parameter values Takeuchi et al. (2009). Furthermore, we compare chimera states in these mean-field coupled SL ensembles with those arising in SL oscillators subject to a nonlinear global coupling that conserves a harmonic oscillation of the ensemble mean. This coupling scheme was introduced to describe experiments on Si electrodissolution Schmidt et al. (2014); Schönleber et al. (2014). In our context, it gives additional information of how a global constraint might further impact the interplay between high-dimensional incoherent motion and collective interaction patterns.

The paper is organized as follows. In the next section, we introduce our two models as well as the particular types of chimera states that are subject to our investigations. In section III, the notation used throughout the paper is introduced, Lyapunov analysis as performed on our data reviewed, and measures used to characterize the dynamics, in particular the Lyapunov dimension $D_{\text{L}}$ and the inverse participation ratio (IPR) are defined. In the results and discussion section IV, we start by investigating the Lyapunov spectra of the two types of chimera states for $N=16$ oscillators in detail, discuss then how the spectra change when $N$ is stepwise doubled up to $N=256$, and end with a discussion of how the Lyapunov dimension and the IPR depend on $N$. The latter allows us in particular to draw conclusion on the existence and on some properties of collective Lyapunov modes. The paper is completed with the conclusion section V.

II Dynamical systems

There is a vast number of different chimera states that have been reported in literature in recent years Panaggio and Abrams (2015); Kemeth et al. (2016); Schöll (2016); Omel’chenko (2018). In this article, we focus on two kinds, also called type I and type II chimeras Schmidt and Krischer (2015), which can be observed in systems with long-range interactions Sethia and Sen (2014); Schmidt et al. (2014). The former, the type I chimera, appears in Stuart-Landau ensembles with linear global coupling. The time evolution of the systems is governed by equations of the form

[TABLE]

wherein $W_{k}$ denotes the complex amplitude of oscillator $k$, and $k=1,\dots,N$, with $N$ being the total number of oscillators. Parameters are the real shear, $c_{2}$, and the complex coupling constant, $\kappa$, which we set to $c_{2}=2$ and $\kappa=0.7\,(1-i)$ in the following. Hereby, $i$ denotes the imaginary unit. Given this set of parameters, we investigate ensemble sizes $N=2^{l}$ for $4\leq l\leq 8$, and find the type I chimera state to be a stable attractor for $N\geq 16$. The oscillator dynamics in the complex plane, as well as the corresponding time series of the real parts are shown in Fig. 1 (a) and (c), respectively, for $N=256$. There, one observes a coherent group of oscillators, depicted in blue, with a larger absolute value of $W_{k}$, and an incoherent group composed of oscillators with smaller amplitudes, shown in red. The thick dots highlight a snapshot of the dynamics. Note that the motion of the incoherent group, when viewed in the complex plane, lies on a string-like object undergoing stretching and folding. The structure is similar to a Birkhoff-Shaw attractor Thompson and Stewart (1984), and resembles the chaotic motion found in mean-field coupled SL ensembles Nakagawa and Kuramoto (1994); Takeuchi et al. (2009).

In contrast to the type I chimeras, the type II chimera state has only been observed in systems with nonlinear global coupling Schmidt and Krischer (2015). In particular, it appears in Stuart-Landau ensembles of the form

[TABLE]

with parameters $c_{2}=-0.66$, $\nu=0.1$ and the initial absolute value of the mean amplitude,

[TABLE]

Exemplary simulation data are depicted in Fig. 1 (b) and (d) for $N=256$. The coloring scheme is identical to that of the type I chimera state. Note that the type II chimera state exhibits an additional frequency component in the oscillator dynamics and that there is no clear amplitude separation between the coherent and incoherent oscillators, as opposed to the type I dynamics. Furthermore, qualitative differences between both types of chimera states can be captured in terms of order parameters Kemeth et al. (2016): While the type II chimera states possess an order parameter with oscillatory behavior, leading to categorization as ’breathing chimeras’, the order parameter of the type I chimera state fluctuates irregularly around a constant value, as characteristic for the class of ’turbulent chimeras’. In the following, we investigate the ensemble dynamics of the two cases more carefully using Lyapunov analysis, with a particular emphasis on the apparently chaotic motion of the oscillators in the incoherent groups.

III Background and methods

III.1 Lyapunov analysis

Lyapunov analysis provides quantitative methods to investigate the degree of chaoticity of dynamical systems Pikovsky and Politi (2016). To introduce concepts and notation, we consider a sufficiently well-behaved dynamical system, whose time evolution is described by a time-dependent vector $\bm{x}=\bm{x}\!\left(t\right)\in\mathbb{R}^{2N}$, and governed by the autonomous ordinary differential equation

[TABLE]

Eqs. (1) and (2) can be cast in this form by introducing the real-valued oscillator coordinates $\left\{a_{k},b_{k}\right\}_{k=1}^{N}$ through

[TABLE]

wherein $W_{k}^{*}$ denotes the complex conjugate of $W_{k}$. Organizing all coordinates in a $2N$-dimensional state vector then yields

[TABLE]

Furthermore, we specify a reference trajectory $\hat{\bm{x}}\!\left(t\right)$, which at any point in time has to satisfy Eq. (3), and investigate the growth or decay of small perturbations around $\hat{\bm{x}}\!\left(t\right)$. The time evolution of infinitesimal perturbations $\delta\bm{x}=\delta\bm{x}\!\left(t\right)\in T_{\hat{\bm{x}}\left(t\right)}\mathbb{R}^{2N}$ is governed by the tangent-space dynamics

[TABLE]

wherein $\left.\frac{\partial\bm{f}}{\partial\bm{x}}\right\rvert_{\hat{\bm{x}}\left(t\right)}$ is the Jacobian matrix of $\bm{f}\!\left(\bm{x}\right)$, evaluated at $\hat{\bm{x}}\!\left(t\right)$, and $T_{\hat{\bm{x}}\left(t\right)}\mathbb{R}^{2N}$ denotes the tangent space of the dynamical system in the vicinity of $\hat{\bm{x}}\!\left(t\right)$. For our purposes, we have $T_{\hat{\bm{x}}\left(t\right)}\mathbb{R}^{2N}\cong\mathbb{R}^{2N}$.

Given a suitable initial condition $\delta\bm{x}\!\left(t_{0}\right)=\delta\bm{x}_{0}\in\mathbb{R}^{2N}$, Eq. (4) possesses the analytic solution

[TABLE]

wherein $\Psi_{\hat{\bm{x}}}\!\left(t,t_{0}\right)$ denotes the fundamental matrix of Eq. (4), encoding the time evolution of perturbation vectors from $T_{\hat{\bm{x}}\left(t_{0}\right)}\mathbb{R}^{2N}$ at time $t_{0}$ to $T_{\hat{\bm{x}}\left(t\right)}\mathbb{R}^{2N}$ at time $t$. Using Eq. (5), as well as the standard scalar product in $\mathbb{R}^{m}$, the squared amplitude of some perturbation vector $\delta\bm{x}\!\left(t\right)$ can thus be expressed as

[TABLE]

Oseledets’ theorem, also known as the multiplicative ergodic theorem Oseledets (1968), states that under rather general conditions there exist symmetric matrices

[TABLE]

with real positive eigenvalues $\mu_{1}>\mu_{2}>\ldots>\mu_{r}$, $r\leq 2N$, which thus characterize the long-term expansion and contraction rates of perturbation vectors. Note here that symmetries of the system as well as of the dynamical pattern can give rise to degeneracies in the spectrum, which we will discuss later on in more detail. Denoting with $g_{k}$ the degree of degeneracy of eigenvalue $\mu_{k}$, the degeneracies satisfy $\sum_{k=1}^{r}g_{k}=2N$. Thereby, the degeneracies $g_{k}$ also reflect the dimensionality of the eigenspaces $U_{\hat{\bm{x}}\left(t_{0}\right),k}^{(\pm)}$ of $M_{\hat{\bm{x}}\left(t_{0}\right)}^{(\pm)}$, associated with eigenvalue $\mu_{k}$. Furthermore, the eigenvalues $\mu_{k}$ and degrees of degeneracy $g_{k}$ coincide for $M_{\hat{\bm{x}}\left(t_{0}\right)}^{(+)}$ and $M_{\hat{\bm{x}}\left(t_{0}\right)}^{(-)}$, and are independent of the particular reference trajectory for ergodic dynamics. The Lyapunov exponents $\lambda_{1}>\lambda_{2}>\ldots>\lambda_{r}$ are then defined through $\lambda_{k}=\ln\mu_{k}$ and satisfy

[TABLE]

for $\delta\bm{x}\in\Omega_{\hat{\bm{x}}\left(t_{0}\right),l}$. Therein, the subspaces $\Omega_{\hat{\bm{x}}\left(t_{0}\right)}^{(k)}\subset T_{\hat{\bm{x}}\left(t_{0}\right)}\mathbb{R}^{2N}$ are obtained from suitable intersections of eigenspaces of $M_{\hat{\bm{x}}\left(t_{0}\right)}^{(+)}$ and $M_{\hat{\bm{x}}\left(t_{0}\right)}^{(-)}$. Employing the above notation, we define

[TABLE]

with $\oplus$ denoting the direct sum of vector spaces, and obtain

[TABLE]

The set of subspaces $\left\{\Omega_{\hat{\bm{x}}\left(t_{0}\right),k}:1\leq k\leq r\right\}$ is called Oseledets’ splitting Oseledets (1968); Ginelli et al. (2007) and provides a non-orthogonal decomposition of the tangent space according to different expansion rates of infinitesimal perturbations, i. e.

[TABLE]

The Oseledet subspaces satisfy $\dim\left(\Omega_{\hat{\bm{x}}\left(t_{0}\right),k}\right)=g_{k}$. In contrast to the eigenspaces $U_{\hat{\bm{x}}\left(t_{0}\right),l}^{(\pm)}$ of $M_{\hat{\bm{x}}\left(t_{0}\right)}^{(\pm)}$, Oseledets’ splitting is norm-independent, invariant under time inversion, and depends on the current system state in a way that is covariant with respect to the dynamical flow, i. e.

[TABLE]

The spanning vectors of Oseledets’ splitting are called covariant Lyapunov vectors (CLVs) Ginelli et al. (2007); Kuptsov and Parlitz (2012); Ginelli et al. (2013), and by virtue of Eq. (7) indicate the local orientation of stable and unstable manifolds in phase space. The spanning vectors of the orthogonal eigenspaces of $M_{\hat{\bm{x}}\left(t_{0}\right)}^{(\pm)}$ are known as forward (FOLVs, $(+)$) and backward orthogonal Lyapunov vectors (BOLVs, $(-)$). Despite the superior dynamical properties of CLVs compared to FOLVs and BOLVs, the efficient computation of CLVs has become possible only recently due to algorithms by Ginelli et al. Ginelli et al. (2007), Wolfe and Samelson Wolfe and Samelson (2007), and later on Kuptsov and Parlitz Kuptsov and Parlitz (2012), who improved on the method by Wolfe and Samelson. CLVs have since attracted a large amount of scientific interest and proved to be a fruitful source of insight, especially into phase-space structures of high-dimensional dynamical systems. In particular, CLVs have been used to study e. g. dynamics of rigid disk systems Bosetti and Posch (2010); Truant and Morriss (2014), chaotic motion in spatially extended systems Yang and Radons (2013); Xu and Paul (2016), stability properties of geophysical models Vannitsem and Lucarini (2016) and collective chaos in systems of coupled oscillators Takeuchi et al. (2009, 2011); Takeuchi and Chaté (2013).

III.2 Lyapunov dimension and inverse participation ratio

The dimensionality of a chaotic attractor can be estimated from the Lyapunov dimension

[TABLE]

with $L\leq 2N$ denoting the largest integer for which $\sum_{l=1}^{L}\lambda_{l}>0$. Note here that for computing $D_{\text{L}}$, we have to take account for degeneracies in the Lyapunov spectrum explicitly. The summation in Eq. (8) therefore runs over an extended index $l\in\left\{1,\ldots,2N\right\}$, which enumerates all of the $2N$ potentially degenerate LEs in a way that satisfies $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{2N}$. Under generic circumstances, the Kaplan-Yorke conjecture then states that $D_{\text{L}}$ is an upper bound of the information dimension $D_{\text{I}}$ of the dynamical pattern Kaplan and Yorke (1979), which is obtained as a special case of the Renyi-$q$ dimension for $q=1$ Rényi (1970); Grassberger (1983), and is closely related to the information production in the underlying system Ott (1997). A linear scaling of $D_{\text{L}}$ with the system size, i. e. in our case the number of oscillators $N$, is commonly used to demonstrate extensivity of chaotic dynamics.

As mentioned above, oscillatory systems with global coupling schemes might also possess collective modes arising from strong correlations Matthews et al. (1991); Hakim and Rappel (1992); Nakagawa and Kuramoto (1993), which do not scale linearly with system size. Such collective modes can be identified investigating the localization or delocalization properties of CLVs associated with particular LEs. Using oscillator coordinates and the above notation, an arbitrary perturbation can be expressed as

[TABLE]

wherein $\delta a_{k}$ and $\delta b_{k}$ denote the relative perturbation amplitudes affecting the real and imaginary part of the oscillator state $W_{k}$, respectively. As demonstrated by Takeuchi et al. Takeuchi et al. (2009), the IPR defined by Mirlin et al. Mirlin (2000) is a suitable measure for vector localization. For perturbations in the form of Eq. (9), the IPR can be written as

[TABLE]

Therein, we assume vector normalization according to $\sum_{k=1}^{N}\left(\delta a_{k}^{2}+\delta b_{k}^{2}\right)=1$. By definition, the IPR then takes values between $N^{-1}$ and $1$. Large values are obtained if the vector under study possesses a small number of large-amplitude components, indicating localization of the perturbation mode. Smaller values are obtained if all terms in the summation are of similar magnitude so that a greater number of oscillators is affected by the perturbation. A small IPR value is thus indicative of a collective Lyapunov mode.

III.3 Numerical methods and simulation details

For the integration of the system dynamics, we use the variable-step Dormand-Prince method Dormand and Prince (1980) implemented in the explicit Matlab integration function ode45 Shampine and Reichelt (1997) with a maximum time step of $5\cdot 10^{-3}$ for Eq. (1) and $1.5\cdot 10^{-2}$ for Eq. (2). During simulation of the type I chimera states, we allowed an initialization period of at least $10000N$ periods of the dynamics for the system state to settle down to an attractor. In the case of type II chimera states, we observe extremely long oscillatory transients, during which the number of synchronized oscillators increases with time. The duration of the transients is found to scale exponentially with system size and is of the order of $10^{6}$ periods of the dynamics for $N=256$. For a give system size, we find that the time between successive oscillators joining the synchronized cluster grows exponentially with increasing size of the cluster. The observed dynamics might be related to super-transients as reported in Wolfrum et al. (2011), but a more detailed analysis of the transients is beyond the scope of this paper. In order to minimize the influence of transient dynamics on our studies, we chose the initialization period so long as to observe no more changes in the number of synchronized oscillators for at least $10000N$ periods. Based on an analysis of Fourier spectra, we assume average period lengths of the dynamics of approximately $4.2$ time units for type I and $12.9$ for type II chimera states, independent of the system size $N$. For the computation of LEs and BOLVs, we employ the standard method by Eckman et al. Eckmann and Ruelle (1985), which is based on repeatedly performing $QR$-decompositions and averaging over the logarithmic diagonal entries of the $R$ matrices for a sufficiently long time. We choose to perform 20 $QR$-decompositions per period of the dynamics. Expecting an exponential convergence of the $Q$-matrices containing the BOLVs, we admit a transient period of at least $1000N$ periods before recording the LEs and BOLVs. The CLVs are obtained from the BOLVs and the respective $R$ matrices by applying the dynamic algorithm by Ginelli et al. Ginelli et al. (2007). In backward-time direction, we admit an initialization period of at least $100N$ periods before recording the CLVs. The final estimate of the LEs is obtained from an average over at least $2.5\cdot 10^{4}$ periods of the dynamics. The CLVs are saved after each $QR$-step (resulting in 20 samples per period) over a time interval of $2.5\cdot 10^{4}$ periods. Based on the statistics of the short-time estimates, we expect an accuracy of the exponents of at least $10^{-3}$.

IV Results and Discussion

IV.1 Lyapunov spectra for $N=16$ oscillators

Fig. 2 shows Lyapunov spectra of the type I and type II chimera states for a system size of $N=16$ oscillators. Simplifying further discussions, we plot the full spectra, consisting of $2N=32$ exponents, without removing degenerate exponents. We furthermore follow the suggestion of Pikovsky and Politi (2016) and introduce the rescaled index

[TABLE]

which enables us to show spectra for different ensemble sizes over the same abscissa range. The simulation data, used to generate the spectra, is not depicted here, but resembles the data in Fig. 1 closely. For the type I dynamics, we find a group of $N_{\text{sync}}=10$ synchronized oscillators, which coexists with $N_{\text{inc}}=6$ oscillators performing incoherent motion. For the type II dynamics, the fraction of synchronized oscillators is higher with $N_{\text{sync}}=13$ and $N_{\text{inc}}=3$. The clustering of the oscillators into synchronized and incoherent groups induces a similar grouping pattern also for the LEs: These can be assigned to four main groups – P, N and $\text{CI}_{1\text{,}2}$ –, which are separated by a small number of singleton exponents – Z and $\text{S}_{1\text{,}2}$.

Starting with P, the group of positive exponents, we recognize that both types of chimera possess $N_{\text{inc}}-1$ positive exponent, revealing their hyperchaotic character. In addition to P, we find a group N of negative exponents, with the same number of exponents, cf. Fig. 2. The facts that P and N are identical in size, and that the number of exponents appears to relate to the number of incoherent oscillators in the ensemble, suggests that P- and N-exponents correspond to stable and unstable directions in tangent space, which affect the group of incoherent oscillators predominantly. The shapes of the corresponding CLVs support this finding. In particular, all CLVs associated with the P- and N-exponents perturb all synchronized oscillators in the same way, without disturbing the clustering pattern.

Besides groups P and N, we find a group Z of very small exponents with values of the order of $10^{-4}$ or smaller. We observe two such exponents for the type I chimera state and three for the type II dynamics. These numbers are consistent with the number of zero-exponents to be expected from symmetry considerations of Eqs. (1) and (2). The governing laws of both systems are independent of time, so that for any solution $\hat{\bm{W}}\!\left(t\right)$ of Eqs. (1) or (2) the time-shifted solution $\hat{\bm{W}}\!\left(t+\delta t\right)$ provides a valid trajectory, as well. Perturbations along the dynamical flow therefore neither grow nor decay in time, resulting in a neutrally stable direction of perturbation and a zero LE. Trivially, the invariance is preserved also when changing to the real-valued coordinates. With Eq. (3) it is easy to see that the corresponding CLV is

[TABLE]

which, by chain rule, satisfies Eq. (4). Similarly, Eqs. (1) and (2) are invariant with respect to phase shifts in the complex plane, i. e. for any angle $\phi\in\mathbb{R}$, $e^{i\phi}\,\hat{\bm{W}}\!\left(t\right)$ is a valid solution if the same is true for $\hat{\bm{W}}\!\left(t\right)$. As a result, another zero-exponent is associated with the covariant Lyapunov vector

[TABLE]

corresponding to an infinitesimal rotation of all oscillators in the complex plane.

For the type II chimera state in Eq. (2), we have an additional zero exponent, which arises from the conservation law

[TABLE]

This can be seen from the fact that

[TABLE]

wherein $\left\langle W\right\rangle=\sum_{k=1}^{N}W_{k}/N$ denotes the complex amplitude of the mean field. According to Eq. (12), $\left\langle W\right\rangle$ performs a harmonic oscillation in time, and thus conserves the mean-field amplitude. The conservation law (11) induces a splitting of phase space into invariant manifolds associated with different values of $\left\lvert\left\langle W\right\rangle\right\rvert=\eta$, and thus yields another neutrally stable direction in tangent space. The CLV associated with this direction, however, appears to possess a non-trivial structure so that an analytical expression has not been identified, yet.

Two further groups of exponents, $\text{CI}_{1,2}$, are made up of only two distinct LEs, possessing a $\left(N_{\text{sync}}-1\right)$-fold degeneracy, each. Clearly, these degeneracies originate from the synchronization pattern of the chimera states. To see this, we have to take a closer look on the tangent space dynamics in presence of synchronized oscillators. Extending the notation of Eq. (3), we write

[TABLE]

with $\bm{a}=\left(a_{1},\ldots,a_{N}\right)^{T}$ and $\bm{b}=\left(b_{1},\ldots,b_{N}\right)^{T}$. The Jacobian matrix of the overall system can then be written in block-matrix form as

[TABLE]

Denoting the $k^{\text{th}}$ coordinate of $\bm{f}^{(a)}$ with $f_{k}^{(a)}$, the global coupling scheme in Eqs. (1) and (2) allows us to separate the terms according to

[TABLE]

with $c^{(a)}$ and $g^{(a)}$ being independent on the oscillator indices $k$ and $j$, respectively. A similar argument is valid for $\bm{f}^{(b)}$. Sorting then the oscillators in a way that synchronized oscillators $W_{k}$ obtain indices $k=1,\ldots,N_{\text{sync}}$, we can conclude from Eq. (14) that each of the sub-Jacobians in Eq. (13) can again be written in block-matrix form as

[TABLE]

wherein for $p,q\in\left\{a,b\right\}$,

[TABLE]

and $\bm{1},E\in\mathbb{R}^{N_{\text{sync}}\times N_{\text{sync}}}$ denote the unit matrix and the matrix with entries all equal to unity. Furthermore, the entries of $R_{pq}\in\mathbb{R}^{N_{\text{inc}}\times N_{\text{sync}}}$ are row-wise constant. From this structure, we can deduce that for perturbations affecting only the relative position of synchronized oscillators, i. e.

[TABLE]

with $\sum_{k=1}^{N_{\text{sync}}}\delta a_{k}=\sum_{k=1}^{N_{\text{sync}}}\delta b_{k}=0$, we have

[TABLE]

as well as

[TABLE]

for $k=N_{\text{sync}}+1,\ldots,N$. Perturbations of type (15) thus form an invariant subspace under time evolution in tangent space, the dimension of which is $2N_{\text{sync}}-2$. Reducing tangent space dynamics to only this subspace, it is possible to obtain a two-dimensional system of linear equations, whose temporal asymptotics exactly reproduce the values of the two degenerated LEs, as obtained from the full system dynamics. Similar sets of exponents have also be found by Ku et al. Ku et al. (2015) in the context of pure cluster states in oscillatory systems with linear global coupling, and were termed cluster integrity exponents. Due to the very similar role of the degenerate exponents in our context, we adopt this nomenclature here, leading to the group annotations $\text{CI}_{1}$ and $\text{CI}_{2}$.

What remains in both spectra is a small number of singleton exponents – $\text{S}_{1}$ and $\text{S}_{2}$ in the type I spectrum and only $\text{S}_{1}$ for type II –, the role of which will become clear when investigating the localization properties of the associated covariant Lyapunov vectors.

IV.2 Lyapunov spectra for larger systems

Fig. 3 shows Lyapunov spectra for type I and type II chimera states in Eqs. (1) and (2). Again, the spectra are plotted against the rescaled index and include the degenerate exponents. Apart from small fluctuations of the order of two to four oscillators, the number $N_{\text{inc}}$ of incoherent oscillators scales roughly linearly with ensemble size, so that the structure of the spectra is similar to those in Fig. 2. In particular, we find that in most cases the number of P- and N-exponents equals $N_{\text{inc}}-1$ each, and thus scales extensively with system size. Only for the type I spectra with $N=128$ and $N=256$ oscillators, there are $N_{\text{inc}}-2$ positive exponents and one additional singleton exponent very close to zero. This, however, is likely to be an artifact. On the one hand, especially for high-dimensional systems, it may take a long time for the smallest exponents to converge properly, since these correspond naturally to the slowest directions of perturbation. On the other hand, the $QR$-based method for computing LEs is known to be vulnerable to numerical errors in computing close-to-zero exponents and associated BOLVs Dieci et al. (1997); Dieci and van Vleck (2005); Pikovsky and Politi (2016).

Similar to P- and N-exponents, the groups of degenerate exponents scale extensively with system size, as well. In agreement with theory, we find a degeneracy of degree $N_{\text{sync}}-1$ for each of the exponents.

While the values $\lambda_{l}$ of the degenerate exponents do not appear to follow an obvious trend as a function of system size, the P- and N-exponents exhibit a pronounced size dependence with a significant trend. To give a better overview on the size dependence of the non-degenerate P- and N-exponents, the smaller sub-figures in Fig. 3 depict magnified views onto the relevant regions. For the type I dynamics, we observe a decreasing trend for the main part of the positive exponents while some of the positive exponents close to the zero group appear to follow an increasing trend. The bulk of P-group thus seems to approach a constant value, which is consistent with former results by Takeuchi et al. Takeuchi and Chaté (2013). However, for the most positive exponents, the rate of decrease is slower compared to bulk and furthermore appears to slow down with increasing $N$. Similar to the observations of Takeuchi et al., this might indicate a sub-partitioning of the positive exponents into an extensively scaling and a sub-extensive group.

In contrast to the type I case, we observe an increasing trend for the main part of the positive exponents within the type II spectra. Only some of the smallest positive exponents do not follow the increase. In particular, the smallest positive exponent seems to approach zero, so that for $N=256$ oscillators four instead of three exponents are very close to zero. From a numerical perspective, the values of the Z-group exponents of the type II dynamics fluctuate stronger compared to the those of the type I spectra. This might indicate difficulties with the numerical integration method. Still, it might also give an indication that in the limit of large system sizes some of the CLVs come very close to those of the zero subspace and should thus be explored further.

The N-group of negative exponents appears to flatten with increasing system size for both types of chimera states, i. e. more negative exponents increase whereas less negative exponents decrease. For the type I spectra, this results in an increasingly tangent-like shape of the spectral curve. A similar pattern is found also within the type II spectra. A difference, however, occurs in the behavior of the most negative exponents between $\tilde{l}=0.95$ and $\tilde{l}=1$. Within this range, the type II spectra deviate significantly from the tangent-like shape. The presence of an additional turning point in the spectral curve around $\tilde{l}=0.975$ results in an step-like shape of the N-group spectrum and might suggest that there exist additional substructures in the spectrum, which are not resolved by the grouping structure proposed in this paper.

IV.3 Lyapunov dimension and extensivity of type I and type II chimera states

Given the Lyapunov spectra, it is straight forward to compute also the Lyapunov dimension as stated in Eq. (8). The computed values are shown in Fig. 4 (a) and reveal a monotonic increase of $D_{\text{L}}$ with system size. Interestingly, the trend of the type I chimera state (at least for system sizes below $N=256$) is not perfectly linear, but slightly slower. Going to larger system sizes, however, the deviation from linear growth decreases, suggesting that the attractor dimension of the chimera state behaves extensively in the limit $N\to\infty$. Similarly, the Lyapunov dimension appears to grow slightly faster than linear for type II chimera states, but growth slows down when considering larger systems.

The extensive nature of the attractor dimension is further strengthened by another consideration. Due to identical synchronization of $N_{\text{sync}}=N-N_{\text{inc}}$ oscillators, the system state $\hat{\bm{x}}\!\left(t\right)$ is necessarily confined to a $2(N_{\text{inc}}+1)$-dimensional subspace of the $2N$-dimensional phase space. From Eqs. (1) and (2), it is easy to see that oscillators remain synchronized for all time, once they move in synchrony, so that this subspace is invariant under the dynamical flow. Although simulations achieve synchronization only within numerical accuracy, it thus seems natural to compute the Lyapunov dimension in a way that excludes directions orthogonal to the synchronization manifold. In section IV.1, we have seen that these directions correspond to the cluster integrity exponents $\text{CI}_{1}$ and $\text{CI}_{2}$. Therefore, removing the $2(N_{\text{sync}}-1)$ degenerate exponents before computing $D_{\text{L}}$ results in an effective value of the Lyapunov dimension, $D_{\text{L}}^{\text{(eff)}}$, which respects the synchronization pattern. This procedure is equivalent to considering an $(N_{\text{inc}}+1)$-oscillator system, in which the synchronized oscillators are replaced by a single representative and weighted accordingly before taking averages. When plotting this effective dimension as a function of the reduced system size $(N_{\text{inc}}+1)$, finite-size effects are minimized and the deviations from linear growth are marginal, as shown in Fig. 4 (b). To quantify the deviation, we presumed a power-law dependence of the form

[TABLE]

and obtained power-law exponents $\beta=0.96$ for the type I chimera and $\beta=0.94$ for the type II chimera from a log-log linear fit to the data. Within numerical accuracy, both exponents are well consistent with an extensive scaling of the effective Lyapunov dimension, and thus with extensively chaotic motion of the incoherent oscillators.

IV.4 Localization of covariant Lyapunov vectors

As a last part of the Lyapunov analysis of type I and type II chimera states, we discuss the localization properties of the CLVs in terms of the IPR as stated in Eq. (10). Generically, due to time dependence of the CLVs themselves, we find that the IPR value for a given vector fluctuates strongly as a function of time. This fact has already been noted by Takeuchi et al. Takeuchi et al. (2009), who considered the temporal mean value $\left\langle\text{IPR}\right\rangle_{t}$ in order to overcome this problem. Within our study, we changed this approach slightly and replace the mean with the median of the temporal distribution. Investigating the distribution of IPR values in time, we noticed that the mode of the probability distribution especially for larger ensemble sizes and delocalized modes is typically situated at values that are notably smaller than the temporal mean of the distribution. Increasing the ensemble size, we found that the separation between mode and mean grows monotonically, so that vector delocalization appears to be underestimated systematically by the mean value for larger system sizes. As a reason for this effect, we identified an extremely heavy-tailed shape of the distribution functions, which favors larger values of IPR in a way that introduces bias to the mean. Employing the median of the temporal distribution instead of the mean reduces the dependence on heavy tails and is expected to yield more reliable results. In Fig. 5, we therefore show the temporal median IPR of the CLVs, plotted as a function of the corresponding LEs. Therein, special care had to be taken for the cluster integrity exponents $\text{CI}_{1}$ and $\text{CI}_{\text{2}}$, as well as for the zero-exponents Z. Since all of the CLVs associated with a degenerate exponent share the same invariant subspace in tangent space, and due to linearity of the tangent space dynamics, any linear combination of two or more CLVs is again a valid CLV with the same exponent. The IPR, however, is intrinsically nonlinear and varies strongly when performing linear combinations, so that strictly speaking the IPR can not be properly defined in degenerate subspaces. In practice, of course, well-defined values are obtained, but the values depend on initialization details and thus carry no information on the actual dynamics. For reasons of clarity, we thus excluded the groups $\text{CI}_{1}$ and $\text{CI}_{2}$ from the plots. For the zero-exponents, especially for the type II chimera state, we observed a similar effect. Knowing analytical solutions of the corresponding covariant Lyapunov vectors, however, we chose one vector to be equal to the time derivative of the dynamics $\delta\bm{x}_{\text{ts}}$, another one to equal the phase-shift direction $\delta\bm{x}_{\text{ps}}$, and the remaining one to lie within the zero-subspace, linearly independent of the former ones.

From the localization spectra, we recognize that the grouping of the Lyapunov spectra manifests itself also in the localization properties of the CLVs. Regarding the type I spectra in Fig. 5, we observe a notable decrease of the IPR with increasing system size for the largest few exponents, for close-to-zero exponents and for some of the smallest exponents. These findings are consistent with the results of Takeuchi et al. Takeuchi et al. (2009); Takeuchi and Chaté (2013), who identified five collective Lyapunov modes in similar positions for string-like chaotic states in the same set of equations. Trivially, the zero-exponents, resulting from time- and phase-shift symmetries, are associated with delocalized vectors, yielding two neutrally stable collective Lyapunov modes. Close to zero, however, we observe additional Lyapunov modes that appear to become increasingly delocalized with increasing system size and correspond to both positive and negative near-zero exponents. On the negative side, we observed already in Fig. 2 the close-to-zero singleton exponent $\text{S}_{1}$, which possesses analogues also in spectra for larger system sizes, as visible in Fig. 3 (Type I, P). The exceptional position of this exponent in the spectra might indicate correspondence to some collective Lyapunov mode or to an otherwise special direction in phase space yet to be identified. Similarly, we observe also a small number of positive exponents, for which the median IPR appears to decrease with increasing $N$. However, these close-to-zero modes are hard to separate from the zero-subspace because of both, strong finite-size effects, which were also found by Takeuchi et al. Takeuchi et al. (2009), and the before-mentioned inherent numerical inaccuracies of close-to-zero exponents and associated BOLVs from $QR$-decomposition. We thus leave the study of larger systems, which exceed our current computational facilities, for future research. Nevertheless, the qualitatively similar pattern of Lyapunov spectra and median respectively mean IPRs of the incoherent group of type I chimera states and string-like chaos in ensembles of SL oscillators with linear global coupling reveals that the synchronized group plays only an inferior role for the macroscopic chaotic motion of the ensemble. This could be a hint that the fast collective modes (corresponding to the most stable and unstable CLVs) are only delocalized among the incoherent oscillators while the synchronized group is not affected. Looking at the spectrum of the type II chimera state, it appears flatter than the one for the type I dynamics, and a notable decrease of the median IPR is visible only for zero exponents and the ones close-by, cf. Fig. 5 (Type II). In particular, we note that in the type II spectrum neither the most positive LE nor the singleton exponent $\text{S}_{1}$ appear to be associated with a collective Lyapunov mode, since their median IPR decrease only slightly. This reveals significant differences between the collective properties of type I and type II chimera states. It is tempting to conclude that it is the additional constraint on the dynamics of type II chimeras, the conservation of the mean oscillation (11), that impedes the formation of fast collective motion. However, at the current state of analysis this remains a hypothesis. A starting point for further investigations are type-I-like chimera states which transiently form also in Stuart-Landau ensembles subject to nonlinear global coupling Schmidt and Krischer (2015). If for some parameter values these states can be stabilized or at least pushed to a ’super-transient regime’, a Lyapunov analysis should reveal whether also then the collective modes are suppressed, or whether the string-like chaos of the incoherent group is a signature of its collective dynamics.

V Conclusion

In the present paper, we have investigated Lyapunov spectra and localization properties of the associated covariant Lyapunov vectors (CLVs) for two examples of chimera states in systems of globally-coupled amplitude oscillators. We found that for both type I and type II chimera states the Lyapunov exponents can be subdivided into four main clusters, containing extensively many LEs, and a small number of singleton and zero exponents. Two of the main clusters, $\text{CI}_{1}$ and $\text{CI}_{2}$, are found to consist of degenerate exponents, which are a result of the synchronization pattern of the chimera state. The corresponding perturbations affect the synchronized oscillators only, leading to the name cluster integrity exponents. A further group contains the positive exponents and indicates the hyperchaotic nature of the chimera states. For both type I and type II chimeras, the spectra exhibit a pronounced dependence on system size. While, for the type I dynamics, most of the positive LEs decrease with increasing system size towards a limiting value, the main part of the positive exponents grows to a limiting value in the case of type II dynamics. The origin of this difference in behavior has to be addressed in future research.

Notwithstanding the different trend of the positive LEs of our two types of chimera states, the Lyapunov dimension of both dynamics grows extensively with system size. Next to the standard Lyapunov dimension, computed from the whole spectrum of LEs, we defined an effective Lyapunov dimension by neglecting the cluster integrity exponents, which measures the dimensionality only of the effective dynamics and takes account for the synchronization pattern. The size dependence of the effective Lyapunov dimension can be fitted with a power law, and yields a power-law exponent very close to unity for both type I and type II dynamics.

An analysis of the localization properties of covariant Lyapunov vectors for type I and type II chimera states revealed further differences between the dynamics. While we found evidence for at least four collective Lyapunov modes being present in the type I dynamics, two of them being associated with large positive and negative LEs, the data suggests that no more than four weak collective modes with LEs zero or close to zero exist in the type II chimera state. Future research should set in at this point and provide additional data for larger system sizes to exclude finite size effects from the localization spectra and to determine the exact number of collective Lyapunov modes for both types of chimera state. Here, it would be of particular interest whether the collective modes form a sub-extensive group that grows approximately as $\mathcal{O}\!\left(\ln\!\left(N\right)\right)$ as found for string-like chaos Takeuchi and Chaté (2013).

From a theoretical point of view, it can be expected that constraints on the system dynamics have a qualitative impact on tangent-space geometry. It therefore appears promising to investigate the geometric relation between the collective Lypaunov modes, for example in terms of angles or correlations between the vectors. Such a study could reveal information on the effective degrees of freedom within the ensemble, and might thus guide a way towards a statistical description of nonlinear amplitude oscillators subject to global coupling and constraints. To address this topic, it furthermore appears reasonable to examine in which way the collective modes affect the oscillator distribution in the complex plane. Solutions to these problems, however, have to be based on a thorough statistical analysis of the dynamics and thus require the study of oscillatory ensembles of much larger size.

Acknowledgments

The authors thank S. W. Haugland for fruitful discussions. The project was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) project KR1189/18.