Chimerapedia: coherence–incoherence patterns in one, two and three dimensions

Chimera states, or coherence–incoherence patterns in systems of symmetrically coupled identical oscillators, have been the subject of intensive study for the last two decades. In particular it is now known that the continuum limit of phase-coupled oscillators allows an elegant mathematical description of these states based on a nonlinear integro-differential equation known as the Ott–Antonsen equation. However, a systematic study of this equation usually requires a substantial computational effort. In this paper, we consider a special class of nonlocally coupled phase oscillator models where the above analytical approach simplifies significantly, leading to a semi-analytical description of both chimera states and of their linear stability properties. We apply this approach to phase oscillators on a one-dimensional lattice, on a two-dimensional square lattice and on a three-dimensional cubic lattice, all three with periodic boundary conditions. For each of these systems we identify multiple symmetric coherence–incoherence patterns and compute their linear stability properties. In addition, we describe how chimera states in higher-dimensional models are inherited from lower-dimensional models and explain how they can be grouped according to their symmetry properties and global order parameter.


Introduction
Complex networks comprising many rythmically oscillating constituents are often used to model different selforganization phenomena in physics, chemistry and biology [1][2][3][4]. In some cases, for example, in reactiondiffusion systems satisfying certain conditions on local kinetics and diffusivities [5], the interaction between oscillators assumes a particular form, known as nonlocal coupling, where the coupling extends over a range wider than nearest neighbor (i.e. local) coupling and depends on separation. This type of coupling also arises in many neural field models [6], where the coupling strength between spiking neurons may fall off exponentially or even change sign with separation, as in the case of the Mexican hat coupling.
Identical phase oscillators with a nonlocal coupling are capable of exhibiting states of partial synchrony in which a subset of the oscillators oscillates with a common frequency and phase while the rest oscillate incoherently. This unexpected state, first identified in the work of Kuramoto and Battogtokh [7], has since been called a chimera state [8]. Since then similar states have been observed in a variety of other coupled oscillator systems, both in experiments [9][10][11][12][13][14][15][16] and in numerical simulations, see [17][18][19] and references therein. However, a theoretical understanding of these unusual states and in particular of their stability properties remains scarce, particularly for two-and three-dimensional oscillator arrays. In the present work we study a particular realization of systems exhibiting partial synchrony where chimera states can be computed exactly, albeit semi-analytically, together with their stability properties, and the results compared with the results of direct numerical simulations of large numbers of oscillators [17]. In particular, we show that two-and three-dimensional oscillator arrays can be analyzed using techniques used earlier for one-and two-dimensional systems [20,21], resulting in a substantial advance in our understanding of higher-dimensional chimera states. Moreover, a systematic exploitation of our techniques enables us to uncover a number of new chimera states in higher-dimensional systems, and to group them into families distinguished by their symmetry properties and values of their global order parameter. These results are obtained using a continuum analysis valid in the limit of infinitely many oscillators [22][23][24][25][26].
We consider three systems of identical nonlocally coupled phase oscillators of increasing complexity, oscillators on a ring, on a two-dimensional torus and on a three-dimensional torus. In the first of these cases we study the dynamics of a ring of phase oscillators Y where αä(−π/2, π/2) is the phase lag parameter and is a 1D cosine kernel determining a specific index-dependent coupling between oscillators. The model(1)- (2) was first suggested and analyzed in [8], and then studied in more detail in [27]. However, as will be shown below, its understanding is still far from complete.
If instead of the one-dimensional oscillator array we consider a square lattice of phase oscillators Y with the 2D cosine kernel Model(3)-(4) was previously studied in [28]. This work reported several types of symmetric chimera states and revealed some relations among them, but did not provide any stability analysis for the reported chimera states. The most complicated model considered in this paper is a cubic lattice of coupled phase oscillators Y  with the 3D cosine kernel cos cos cos , . 6 3 To the best of our knowledge the model(5)- (6) has not been studied and all our results concerning its symmetric chimera states are therefore new. Moreover, these appear to be the first rigorous mathematical results available for three-dimensional chimera states which have up to now only been studied numerically [29][30][31][32][33].
Here and in the following we use the term symmetric chimera to refer to partially coherent states characterized by local order parameter functionsw(x), w(x, y) and w(x, y, z) defined below that are invariant with respect to a nontrivial subgroup of the symmetry group of the corresponding coupling functionG. In all three cases we explore the dynamics of our system in the two-dimensional parameter space (α, A).
The choice of the cosine kernels(2), (4) and (6) is motivated by the fact that this is one of the rare cases where a stability analysis of the observed chimera states becomes feasible, see section 2. Moreover, varying the kernel parameterA we can sweep over several qualitatively different coupling topologies. For example, when A=0 the kernel(2) determines a global (all-to-all) coupling in model (1). In contrast, when Aä(0,1) the coupling kernel p -¢ ( ( ) ) G k k N 2 describes nonidentical but positive coupling strengths which decay with growing distance between pairs of oscillators. This type of coupling is usually called nonlocal. A different type of nonlocal coupling is obtained for Î ¥ ( ) A 1, . In this case, the short-range coupling between oscillators remains positive, while the long-range coupling between them becomes negative. Finally, in the limit  ¥ A model(1)-(2) is equivalent (after a suitable time rescaling) to the system(1) with the balanced kernel = ( ) G x x cos considered in [21]. Here the terms balanced and non-balanced refer to the sum total coupling strengths The role of the parameterA in the higher-dimensional kernels(4) and(6) is qualitatively the same: varyingA from zero to infinity one passes in succession from global, to positive nonlocal, sign-changing nonlocal and balanced coupling topologies similar to the one-dimensional case. Performing numerical simulations of models(1)-(2), (3)-(4) and(5)-(6) we found eleven different types of symmetric chimera states, see tables 1 and 2, which are naturally grouped together based on the value of the global order parameter and the inheritance principle formulated below. Note that in all snapshots of the chimera states we use the scaled oscillator positions x k =−π+2πk/N, y l =−π+2π l/N and z m = −π + 2πm/N instead of the integer indicesk, l andm. This convention will help us later to compare the reported coherenceincoherence patterns with their continuum limit counterparts.
The primary difference between the chimera states shown in tables 1 and 2 comes from the fact that they have different degrees of global synchrony. More precisely, for every chimera pattern from table 1 the corresponding global order parameter stays close to zero for all time, while for every pattern from table 2 the value of ( ) R t is well-separated from zero and fluctuates around some positive constant level. Notice that according to the theoretical results of section 3 chimera states with vanishing order parameter can be found for both balanced and non-balanced cosine kernels. In contrast, chimera states with positive global order parameter exist in systems with non-balanced kernels only.
In general, every chimera state found in the one-dimensional model(1)-(2) can be trivially extended to become a solution of the two-dimensional model(3)-(4) and of the three-dimensional model(5)- (6). Similarly, every two-dimensional chimera pattern can be lifted into a solution of the three-dimensional model(5)- (6). We call this simple observation the inheritance principle and use it to group the chimera states in different columns of tables 1 and 2. Each of these columns is named after the symmetry group relevant to the topmost pattern in the column (details regarding this classification are presented in section 3). In section 3 we also show that the stability region of every inherited chimera state is a subset of the stability region of the corresponding lower-dimensional parent chimera state. In particular, an inherited chimera state can be Table 1. Chimera states with vanishing global order parameter in models(1)-(2), (3)-(4) and(5)-(6). everywhere unstable (left column in table 1), have the same stability region as the parent chimera (middle and right columns in table 1), or its stability region can be a proper subset of the stability region of the parent chimera state (left and middle columns in table 2).
The paper is organized as follows. In section 2 we introduce the main mathematical tools used to analyze the existence and stability of chimera states in models(1)-(2), (3)-(4) and(5)-(6) with infinitely many oscillators (i.e. in the continuum limit). In section 3 we report our main results, describing them separately for each of the chimera states we have found. In section 4 we briefly discuss the impact of the coupling kernel in models(1), (3) and(5) on the types of observed chimera states. Some concluding remarks are given in section 5. Finally, in the appendix we present a series of mathematical formulas used to optimize the efficiency of our numerical solvers.

Ott-Antonsen equation
If the number N of oscillators in the model(1), (3) or(5) tends to infinity, then the long-term dynamics of its solutions can be described by the integro-differential Ott-Antonsen equation [17] whereu is an unknown complex-valued function, called the local order parameter, andu denotes the complex conjugate ofu. Depending on the dimensionality of the model (1D for model (1), 2D for model(3) and 3D for model(5)) one has to assume Table 2. Chimera states with non-vanishing global order parameter in the models(1)-(2), (3)-(4) and(5)- (6 The symbol in equation (7) denotes respectively the integral operators , , d . 10 In the following we concentrate in greatest detail on the three-dimensional case of equation (7) with the corresponding integral operator(10).

Self-consistency equation for relative equilibria
All stationary chimera states observed in the model(5)-(6) can be identified with relative equilibria of the threedimensional Ott-Antonsen equation (7), i.e. solutions of the form = is a time-independent amplitude satisfying the inequality  | | a 1, and  W Î . The coherent and incoherent regions of this chimera state are determined by the relations = | | a 1 and < | | a 1, respectively. In [17] it was shown that there exists a universal function with the abbreviation m = W a ( ) i e . 15 i Expression (12) implies that the coherent and incoherent regions of a chimera state written in the form(13) are determined as incoh 3 For the cosine kernel(6), the range of the integral operator on the right-hand side of equation (14) is spanned by seven basis functions , , 1, cos , cos , cos , sin , sin , sin , 16 implying that all solutions to equation (14) can be written in the form Moreover, using (6) it is easy to verify that expression(17) yields a solution to equation (14) if and only if the coefficientsw k solve the system Since the solutions to equation (7) are interpreted as local order parameters, see [17], the global order parameter of the relative equilibrium(13) is given by the formula 3 2 a quantity that is independent of time. Moreover, forw determined by expression (17) we always have = | | R w 1 , see equation (18).
Note that the reduction to finite dimension implied by equation (16) is exact. It is this fact that makes the stability computations that follow tractable.

Stability analysis of relative equilibria
Stability of a stationary chimera state in model(5)-(6) is one-to-one related to the stability of the corresponding relative equilibrium (11), which can be analyzed via a linearization of equation (7). This approach has been successfully applied to the study of chimera states in different one-and two-dimensional models of coupled phase oscillators [20,21,27,34]. In this paper, we reformulate it for the three-dimensional Ott-Antonsen equation (5). For this, we substitute the ansatz into the Ott-Antonsen equation (7) and linearize it with respect to the small perturbation v(x, y, z, t), obtaining We skip the details of the derivation of equation (20) and refer the interested reader to [17,20] for similar calculations. These references also show that the spectrum determined by equation (20) consists of two parts: the essential spectrum : , , , c.c. i , ess 2 3 and the point spectrums pt consisting of isolated eigenvalues of finite multiplicity. To find the eigenvalues λäσ pt we use the ansatz , , e , , e , 22 t t which provides a solution to equation (20) if the eigenvalueλ and the components + - Applying the integral operator to both sides of equation (23), we obtain the nonlocal eigenvalue problem for the components of the eigenmode of the local mean field . It follows that the eigenvaluesλ are either real or occur in complex-conjugate pairs. Recall that for the cosine kernel(6), the integral operator has a finite-dimensional range spanned by the functions ψ k (x, y, z), k=1, K, 7, defined in (16). This means that we can solve equation (24) for the perturbation = + - T of the local mean field as the linear combination Substituting this ansatz into equation (24) yields a system of nonlinear equations for the seven pairs of complex coefficientsV k , k = 1, K, 7. Collecting these coefficients into a single vector  Î V 14 , we can rewrite these equations as an equivalent matrix equation where we solve for the eigenvalueλ and the corresponding kernel vector  Î V 14 . The matrixB(λ) has the structure l l l l l l l  for the other harmonic components, m=2, K, 7. Note that according to(13) the matrixB(λ) results from substituting = (| | ) a h w w 2 into expressions(28) and (29) and so only depends on the solution(μ, w) of the selfconsistency equation (14). Importantly, because of the different definitions (28) for all  n 2 and ¹ A 1. The eigenvaluesλ can be found as solutions of the characteristic equation 14 where I n denotes the n×n identity matrix. If all solutions l ¹ 0 to equation (30) lie in the left half-plane, l < Re 0, then the corresponding relative equilibrium(11) is stable. In contrast, if equation (30) has at least one solutionλ=λ * such that * l > Re 0, then the relative equilibrium(11) is unstable.

Results
In this section we report our main results. For every chimera state in tables 1 and 2 we write the corresponding relative equilibrium of the Ott-Antonsen equation (7), analyze its symmetries and study its stability following the method of section 2.3.

Antiphase chimeras in 1D, 2D and 3D
The chimera state shown in the top-left panel of table 1 comprises two equidistant coherent regions of equal size, which are in antiphase to each other. In the following we call this state the antiphase chimera. Self-consistency equation. Antiphase chimeras correspond to solutions of the self-consistency equation (14) given by the formula The minimal dimension needed for chimera states of this type is1D, but they are also inherited in2D and3D cases. Substituting ansatz(31) into (18)- (19) shows that all but one of the equations are automatically satisfied and that the only non-trivial equation reads 1, expression(32) delivers the corresponding values ofΩ andα, see (15). Thus, we obtain an explicit parametric representation of antiphase chimera states. Note that the absence of a constant term in expression(31) analogous tow 1 in(17) implies the global order parameter of every antiphase chimera state vanishes.
Symmetries. The function(31) is invariant under the following discrete symmetry operations: Both these operations are elements of order two (k = 1 , and therefore generate the Klein four-groupD 2 . The symmetry group of antiphase chimeras is therefore D 2 . In this statement we omit the translation and reflection invariance of w with respect to the variables y and z. Stability. In view of expressions (28), (29) and(31) the matrixB(λ) has the following structure (see appendix A.2): The characteristic equation (30) therefore decouples into four independent equations: Recall that antiphase chimeras are 1D geometric patterns that are inherited by the 2D and 3D coupled oscillator models. As a result these states may have two types of instabilities: longitudinal instabilities relevant to the main pattern direction (in the case of expression (31) this is the x-direction), and transverse instabilities appearing for inherited patterns in higher-dimensional models. The longitudinal instabilities are described by equations (33), (34) and (36), while the transverse instabilities appearing in the 2D and 3D cases are determined by equation (35).
Because the blocks of the matrixB(λ) corresponding to the y-instabilities and those corresponding to the z-instabilities coincide, the 2D inherited antiphase chimera and the 3D inherited antiphase chimera are stable or unstable simultaneously. Equation (33) has two types of unstable eigenvalues, complex eigenvalues corresponding to Hopf modes (H) and real eigenvalues corresponding to symmetry-breaking modes (SB). These modes are present for small and large values of the parameter A, respectively, as shown in figure 1 shows the corresponding stability region in the (A, π/2−α) parameter plane. This region is bounded by curves of Hopf (dashed) and symmetry-breaking (dotted) bifurcations and represents the stability region for antiphase chimera states in 1D. In fact, in the present case the complex conjugate eigenvalues emerge from the essential spectrum on the imaginary axis, i.e. from the edge of the essential spectrum, in a manner that is reminiscent of the eigenvalue behavior in the so-called edge bifurcation [35]. For larger values ofA the characteristic equation (33) has instead a real unstable point eigenvalue. This eigenvalue also emerges from the essential spectrum asA increases, although it is subsequently reabsorbed by it. Thus the stability region of antiphase chimeras turns out to consist of two disjoint domains.
In contrast the blocksB nn (λ) with n > 1 can only produce A-independent instabilities. Indeed, expression(15) together with equation (32) imply that the ratio Ω / A does not depend onA. Consequently, if we consider the matrix l¢ ( ) B A nn , we easily verify that it depends onl¢ andp but not onA. In other words, if a blockB nn (λ) with n>1 determines an unstable eigenvalue for someA then this fact remains true for any other ¹ A 0. As a consequence for n>1 it is enough to analyze instabilities produced by blocksB nn (λ) with A=1. In particular we find that equation (34) yields no unstable eigenvalues, although it always determines one zero eigenvalue. Another zero eigenvalue always appears as a solution of equation (36). Moreover, equation (36) has a real unstable eigenvalue for α<α 1 where α 1 ≈1.396, see figure 1(c). This fact accounts for the horizontal line of symmetry-breaking bifurcations in figure 1(a). Figure 1(c) also shows that equation (35) has a real unstable eigenvalue for all αä(0, π/2). This means that antiphase chimeras are always unstable in 2D and 3D and so can only be observed as stable patterns in the 1D case.
3.2. Twisted chimeras (2D) and twisted planes (3D) Self-consistency equation. The twisted plane chimera, figure 2, corresponds to the solution of the self-consistency equation (14) given by the expression It has a 2D counterpart, namely the twisted chimera, but no 1D counterpart. Substituting ansatz (37) into (18)- (19) shows that some of these equations are satisfied automatically, while the remaining ones are all equivalent to a single equation of the form e e e cos . 38  Expression (38) provides an explicit parametric representation of twisted plane chimeras, determiningΩ andα as functions of Î ¥ ( ) p 1 2, , see (15). As for antiphase chimeras, the global order parameter of every twisted plane chimera vanishes.
Symmetries. The function (37) is invariant under the symmetry operations y  w x y   :  , ,  , , ,  : , , , as well as the continuous transformation The first two operations are elements of order two (k = 1 ). Moreover κ 2 and κ 3 do not commute. The symmetry group of twisted plane chimera state is therefore the group Z 2 × O (2). In this statement we omit the translation and reflection symmetries associated with the z-direction.
Stability. In view of (28), (29) and (37)   The characteristic equation (30) therefore decouples into the three independent equations: Notice that equation (40) can be further simplified. Using the column/row interchange property of a determinant we rewrite equation (40) in the form The block symmetry of the latter determinant allows us to represent it as a product of two half-size determinants and to show that every solution of equation (42) solves simultaneously one of the following equations  3.3. Spiral chimera (2D) and spiral rolls (3D) Self-consistency equation. The spiral roll chimera, figure 4, corresponds to a solution of the self-consistency equation (14) given by the expression This chimera state has a 2D counterpart, namely the spiral chimera, but no 1D counterpart. Substituting the ansatz(45) into (18)- (19) shows that some of these equations are again satisfied automatically, while the remaining ones are all equivalent to a single equation of the form Expression(46) provides an explicit parametric representation of spiral roll chimeras, determiningΩ andα as functions of Î ¥ ( ) p 1 2, , see (15). As for antiphase chimeras, the global order parameter of every spiral roll chimera vanishes.
Symmetries. Symmetries of the function(45) were analyzed in [20] in the context of a study of 2D spiral chimeras where it was shown that (45) is invariant under the following symmetry operations: , , e , , . . Therefore the spatial symmetries of the spiral roll chimera constitute the generalized dihedral group of the cyclic groupC 4 , i.e. Dih(C 4 ).   The characteristic equation (30) therefore decouples into five independent equations: In appendix A.4 we explain how one can optimize the computation of the blocks involved in the above equations. Numerical analysis of equations (47)-(51) reveals the stability region of spiral roll chimeras, see figure 5. We find that equations (48), (50) and(51) have simple zero roots for allA andα. Moreover, for small values ofA equation (47) determines a pair of unstable complex conjugate eigenvalues. For fixedα and increasingA these eigenvalues move to the imaginary axis where they are absorbed by the essential spectrum indicating a Hopf bifurcation.
As for equation (36) it can be shown that all instabilities determined by equation (50) are A-independent. Our analysis reveals that for α>α 2 ≈1.23 this equation has a real positive eigenvalue. Hence, the line α=α 2 corresponds to symmetry-breaking bifurcations.
Our extensive search did not reveal any other unstable eigenvalues, including those potentially determined by equations (48) and(49). This indicates that D 3 spiral roll chimeras are stable/unstable for the same parameters (A, α) as their 2D counterparts, i.e. spiral chimeras. Moreover, the stability region from figure 5 apparently extends to arbitrarily large valuesA, a result that agrees with that reported in [34], where stable 2D spiral chimeras were observed in the case of a balanced cosine kernel. 3.4. Classical chimera (1D), coherent stripe (2D) and coherent plane (3D) Self-consistency equation. The coherent plane chimera, figure 6, corresponds to a solution of the self-consistency equation (14) given by the expression cos , where and 0, . 52 0 0 This state has both 1D and 2D counterparts, namely the classical chimera state and the coherent stripe chimera, respectively. Substituting ansatz(52) into (18)- (19) shows that this system is equivalent to the two equations These equations can be understood as follows: we first solve equation (53) forw 0 as a function of p. We then substitute (p, w 0 (p)) into equation (54) and determine the corresponding valuesΩ andα according to expression (15). Note that for the coherent plane chimeras we always have ¹ w 0 0 , implying that for these states the global order parameter is strictly positive, see figure 7.
Symmetries. The function(52) is invariant under the reflection k  -( ·) ( ·) w x w x : , , 1 only and so has Z 2 symmetry.   The characteristic equation (30) therefore decouples into the three independent equations: The stability region of classical chimera states in the 1D case was computed in [27]. It is bounded by curves of fold (solid) and Hopf (dashed) bifurcations, see figure 8. Note that the pair of complex conjugate eigenvalues responsible for the Hopf instability is determined by equation (55). The eigenvalues emerge from the essential spectrum for values ofA greater than one, i.e. for a sign-changing coupling kernel.
The stability regions of coherent plane chimeras (3D) and of coherent stripe chimeras (2D) are identical. They are smaller than the stability region of classical chimeras (1D), because equation (56) has a real unstable eigenvalue for values ofα close toπ/2. This eigenvalue is responsible for a symmetry-breaking bifurcation (dotted curve). Note that λ=0 turns out to be a simple root of both equation (55) and equation (57).

Coherent spot (2D) and coherent tube (3D)
Self-consistency equation. The coherent tube chimera, figures 9(a)-(c), corresponds to a solution of the selfconsistency equation (14) given by the expression This state has a 2D counterpart, namely the coherent spot chimera, but no 1D counterpart. Substituting ansatz(58) into (18)- (19) shows that this system reduces to the simpler set  To solve this system, we first solve equation (59) forw 0 as a function ofp and then use the result to compute the complex quantityμ from equation (60), thereby yielding the corresponding values ofΩ andα, see (15). As for the coherent plane chimeras, the global order parameter of a coherent tube chimera is always strictly positive, see figure 7.
Symmetries. The coherence-incoherence boundary of a coherent tube chimera is determined by the equation = | ( )| w x y z , , 1. According to(58) this yields a cylindrical surface in z with a four-sided squashed circle in the (x, y)-section. It is easy to see that the function(58) has all the discrete symmetries of a square inscribed in the above squashed circle. Hence, all coherent tube chimeras have the symmetry of the dihedral groupD 4 .
Stability. In view of expressions(28), (29) and(58) the matrixB(λ) has the following structure (see appendix A.6): In appendix A.6 we explain how one can optimize the computation of the blocks involved in the above equations. Equations (61) and(63) have simple zeros for all values of the parametersA andα. Moreover, solving equation (61) numerically we found a real positive eigenvalue for smallA and a pair of unstable complex conjugate eigenvalues for largeA. Remarkably, for  A 0.8 the positive real eigenvalue can be found only in a proper subinterval of αä(α fold , π/2), the existence interval for this type of chimera, because at the ends of this subinterval the eigenvalue is absorbed by the essential spectrum. As a result, the stability region of tube chimeras consists of two disconnected parts, see figure 10. In one parameter region the chimera states resemble coherent tubes, figure 9(a)-(c), while in the other region they resemble incoherent tubes, figure 9(d)-(f). Note that in the latter case we can use instead of(58) in order to bring the incoherent tube into the center of the cube p p -[ ] , 3 . The stability region of incoherent tubes lies close to the fold bifurcation curve and is extremely narrow. In contrast, the stability region of coherent tubes lies close to the line α=π/2 and is relatively wide. Notice that all symmetry-breaking bifurcations described by equation (61) are analogous to those known for the coherent spot chimeras in the 2D case [28]. It follows that along with the symmetric coherent plane branch and symmetric coherent tube branch the self-consistency equation (14) also has two additional unstable solution branches corresponding to asymmetric chimera patterns, connecting the symmetry-breaking bifurcation points on the coherent plane branch and the coherent tube branch.
We now turn to equation (62). This equation determines a real positive eigenvalue responsible for another curve of symmetry-breaking bifurcations lying close to the line α=π/2, see figure 10. As a consequence the stability region of coherent tubes in 3D is smaller than the stability region of 2D coherent spots. Since equation (62) has no other unstable roots, the stability regions of 3D incoherent tubes and of 2D incoherent spots are identical.

Coherent ball (3D)
Self-consistency equation. The coherent ball chimera, figure 11, corresponds to a solution of the self-consistency equation (14) given by the expression This chimera state is a purely 3D phenomenon and therefore has neither 2D nor 1D counterparts. Substituting ansatz(64) into (18)- (19) shows that this system reduces to the simpler set To solve this system, we first solve equation (65) forw 0 as a function ofp and then compute the complex quantityμ from equation (66), yielding the corresponding values ofΩ andα, see (15). As for the coherent plane chimeras, the global order parameter of a coherent ball chimera is always strictly positive, see figure 7.
Symmetries. The coherence-incoherence boundary of a coherent ball chimera is determined by the equation = | ( )| w x y z , , 1. According to(64) this yields a six-sided squashed sphere. It is easy to see that the function(64)   In appendix A.7 we explain how one can optimize the computation of the blocks involved in the above equations. Both equations (67) and(68) have simple zeros for all values of the parametersA andα. Moreover, equation (67) has a double real positive eigenvalue for smallA and a double pair of unstable complex conjugate eigenvalues for largeA. Thus the stability region of coherent ball chimeras is bounded by Hopf and symmetrybreaking bifurcation curves as shown in figure 12.
Note that using equations (65)-(66) we can follow the branch of coherent ball chimeras beyond the symmetry-breaking bifurcation. This reveals that the coherent ball chimera transforms continuously into a labyrinthine chimera state, see figure 13. However, the latter coherence-incoherence pattern is unstable and cannot be observed in numerical simulations of equations (5)- (6). On the other hand, such patterns were found in numerical simulations with other coupling kernelsG(x, y, z), see [32].

Other coupling kernels
Multi-harmonic kernels. All of the above results for 3D chimera states can be easily generalized for coupling kernels of the form where  n 2 is an integer and  Î A . Indeed, suppose that (μ, w(x, y, z)) is a solution to the self-consistency equation (14) with the coupling kernel (6). Then (μ, w(nx, ny, nz)) is a solution to the self-consistency equation (14) with the coupling kernel(69) and  n 2. Moreover, comparing expressions (28), (29) written for the solution (μ, w(x, y, z)) and the coupling kernel(6) with their analogs written for the solution (μ, w(nx, ny, nz)) and coupling kernel(69) we find that the resulting stability matricesB(λ) are identical. This leads to the general conclusion: If for some parametersA andα in the model(5)-(6) we observe a chimera state corresponding to a relative equilibrium of the form(13) (for example, one of those shown in tables 1 and 2), then for the same parameters we will find a stable multiple analog of this chimera state in the model(5) with the coupling kernel(69) and  n 2, see figure 14. Note that the sequence of multiple chimera states for multi-harmonic kernels was already observed in 1D and 2D models of coupled phase oscillators [21,34,36]. Here we have shown that there is a fundamental reason for these observations, which we call the chimera multiplication principle.   Mixed harmonic kernels. The Ott-Antonsen equation method described in section 2 can also be used to study chimera states in models(1), (3) and(5) with coupling kernels differing from those considered above. For example, in [21] the authors considered model(1) with the balanced mixed harmonic kernel x  nx  n  x  n  cos  cos  1 ,  .  70 For this model, several new multiple chimera states were found. In terms of the self-consistency equation (14) these correspond to solutions of the form Note that the two-dimensional analogs of these multiple chimeras were also observed [34] in model (3) with the two-dimensional counterpart of the balanced mixed harmonic kernel(70). Therefore according to the chimera inheritance principle we may expect to find three-dimensional analogs of these chimeras too.
Top hat kernels. A large number of numerical results [32,33], including sketches of stability diagrams, were obtained for model (5) with the so-called piecewise constant, or top-hat kernel  p s ps ps where σä(0, 1) is the coupling radius. Along with the chimera patterns from tables 1 and 2 and their multiple analogs, remarkable new states, called knotted and linked chimeras, were found in this model [31,33]. Note, however, that the top-hat kernel has a Fourier series with infinitely many terms, so that the analysis of the corresponding self-consistency equation (14) and especially of the corresponding eigenvalue problem(24) requires a significant computational effort that is beyond the scope of the present work. Importantly, the fact that the top-hat kernel has a Fourier series with infinitely many terms has also another consequence. Similar to the multi-harmonic kernels, every chimera solution(μ, w(x, y, z)) to the selfconsistency equation (14) with the coupling kernel(71) and σä(0, 1) can be transformed into a multiple chimera solution(μ, w(nx, ny, nz)) with  n 2 corresponding to the scaled coupling radiusσ/n. However, the stability properties of the new multiple chimera pattern and of the original chimera pattern may be different, because the eigenvalue problem(24) cannot be reduced to a finite-dimensional equation of the form (30). Hence, the chimera multiplication principle for top-hat kernels works only partially [32,33].

Conclusion
In this paper, we considered basic three-dimensional chimera states observed in the model(5) with the cosine coupling kernel (6). We revealed two important principles, the chimera inheritance principle and the chimera multiplication principle, underlying their appearance. These principles can be easily generalized for higherdimensional situations allowing one to predict the existence and stability of hyper-dimensional chimera states. The simplifications provided by the cosine kernel were used to compute exact stability diagrams of all the reported chimera states and to describe their symmetries. Note that despite the length of the paper our results are far from exhaustive and do not cover all the chimera states that can be found in the model(5)- (6).
First, the self-consistency equation (14) has many symmetric solutions not mentioned in tables 1 and 2. For example, one can use the ansatz cos cos cos , 0, , 72 , , e e e , 0, 73 to construct new interesting relative equilibria as shown in figure 15. However, in numerical simulations of the model(5)-(6) we did not observe any chimera states relevant to these equilibria, and so expect that they are unstable (although this has not been proved rigorously). Second, the numerous symmetry-breaking bifurcations revealed by our analysis indicate that the model(5)-(6) has solutions with less symmetry than those included in tables 1 and 2. This hypothesis is further supported by the observation of asymmetric spot chimeras [28] and asymmetric spiral chimeras [20] in the twodimensional model(3)-(4) as well as by the observation of two pairs of orthogonal spiral rolls in model(5) with the top-hat coupling kernel, see [32]. These states are the result of spontaneous symmetry breaking and must be distinguished from chimera states arising from forced symmetry breaking such as those present in systems with the coupling kernel with distinct A x , A y , A z , B x , B y and B z . Such forced symmetry breaking destroys the highly symmetric chimera states studied in the present work replacing them with states of lower symmetry. Moreover, for nonvanishing coefficientsB x , B y and B z these states typically drift [37]. The lower symmetry of these states makes the corresponding stability calculations more involved, since the linear stability problem may no longer blockdiagonalize. This problem is left for future work.
Third, we have to emphasize that in our work we completely ignored all non-stationary chimera states, for example, quasiperiodic and traveling chimeras, which may potentially appear in the model(5)- (6). The analysis of such states is usually more involved [21,37] and is also beyond the scope of this paper.
Fourth, the solution of the characteristic equation (30) allowed us to determine the stability boundaries of different chimera states but provided no information about the nature of the bifurcations that occur at these boundaries, i.e. whether the reported symmetry-breaking and Hopf bifurcations are subcritical or supercritical and the scaling of the amplitude of the new states generated by these instabilities with the distance from the bifurcation point. Potentially, these questions can be answered via a weakly nonlinear analysis of the Ott-Antonsen equation (7). However, in this case one needs to carry out a non-standard center manifold reduction in a situation where point spectrum eigenvalues emerge from an essential spectrum on the imaginary axis, a situation that leads to significant complications already in the case of globally (all-to-all) coupled oscillators [38][39][40][41].
In view of these complications we have adopted the following numerical approach. For each chimera pattern shown in tables 1 and 2 we carried out a dynamical continuation in four possible directions (increasing or decreasing parameters A and α). For every direction we monitored the time-averaged order parameterR av , see figures 16(a) and (b).
Remarkably, in all cases we observed abrupt jumps in the order parameterR av after the crossing of the corresponding stability boundary. This observation suggests that all the bifurcations reported in this paper are in fact subcritical. More precisely, we found that crossing a symmetry-breaking stability boundary may lead to a collapse of the corresponding chimera state to either a completely coherent state, or to a splay state or another chimera state, which is stable for the new system parameters. In contrast, crossing Hopf instability boundaries generally resulted in abrupt transitions to more complicated time-intermittent regimes, see figure 16(c), suggesting that some of the Hopf bifurcations may be responsible for Type-II intermittent transitions to chaos as described by Pomeau and Manneville in [42]. With other coupling kernels the bifurcations may be supercritical. In particular, supercritical Hopf bifurcations were reported for mixed harmonic [20,21,34], exponential [43] and top-hat [44] kernels, i.e. for kernels containing more than one non-constant Fourier harmonic.
In summary, we have given a systematic (albeit incomplete) overview of three-dimensional chimera states and their relation to the chimera states in lower-dimensional models. The results reported here explain adequately many properties of known symmetric chimera states while motivating further research in the field.

Acknowledgments
The work of OO was supported by the Deutsche Forschungsgemeinschaft under grant OM 99/2-1. The work of EK was supported by the National Science Foundation under grant DMS-1613132.

Appendix. Mathematical tools
A.1. Notation Throughout the appendix we use the following notation for incomplete elliptic integrals while for complete elliptic integrals Note that for n<−1 we use the formula to evaluate the Cauchy principal value of the corresponding singular integral. Moreover, 1 for 0, 0 for 0 denotes the Heaviside function.

A.2. Technical details for antiphase chimeras
The right-hand side of equation (32) can be computed explicitly in terms of elliptic integrals. For this, we first rewrite equation (32) in the form ò ò cos cos d 2 cos cos d . We then divide the integration interval (0, π/2) into two parts (0, x cr ) and (x cr , π/2) where Taking into account the symmetries of antiphase chimeras we find that all off-diagonal blocks in the matrixB(λ) vanish. Moreover, it is easy to verify that   while the corresponding Jacobian reads This allows us to rewrite the self-consistency equation (38)       Recall that in the main text it was shown that equation (40) is equivalent to the two independent equations (43) and (44). Using the above expressions for the blocks B 22 (λ), B 23 (λ), B 25 (λ), B 26 (λ), B 55 (λ) and B 56 (λ) we can simplify the determinants in equations (43) and (44) and show that these equations are equivalent to Thus | | w 2 and w Im 2 are symmetric with respect tox andy, while w Re 2 is antisymmetric with respect tox andy. Using these identities and their symmetries we can show that many off-diagonal blocks of the stability matrixB(λ) vanish, see section 3.3, and that , . 23 32 44 77 To optimize the computation of the right-hand side of equation (46) and of the non-vanishing blocks in the matrixB(λ) we can use the partial integration approach suggested in [20]. All integrals over the square p p Î -( ) [ ] x y , , 2 can be replaced with the double sum of two integrals over smaller squaresD 1 andD 2 . In the squareD 1 we rewrite the integrand using modified polar coordinates In the squareD 2 we rewrite the integrand applying a slightly different coordinate transformation, arcsin cos , arcsin sin , 80 which turns out to have the same Jacobian as in equation (79). Next, we rewrite| | w 2 , x cos 2 and y cos 2 in the new coordinates.     A.6. Technical details for coherent tube chimeras Using the symmetries of the functionw determined by formula(58) we can show that many off-diagonal blocks of the stability matrixB(λ) vanish, see section 3.5, and that l l l l l l l l l l l l To optimize the computation of the right-hand sides of equations (59) and(60) as well as of the nonvanishing blocks in the matrixB(λ) we use the partial integration scheme formulated below. All integrals over the square p p Î -( ) [ ] x y , , 2 can be replaced with the sum of two integrals over smaller squaresD 1 andD 2 . In the squareD 1  In the squareD 2 we rewrite the integrand using a slightly different coordinate transformation, 2arcsin cos , 2arcsin sin , 83 which turns out to have the same Jacobian as in equation (82). Next, we rewrite| | w 2 , x cos 2 and y cos 2 in the new coordinates.
Moreover, for the stability matrixB(λ) we also obtain l l l   To optimize the computation of the triple integrals appearing on the right-hand side of equations (59) and(60) as well as in the definition of the non-vanishing blocks of the matrixB(λ) we use the partial integration scheme formulated below. First, we replace the integration domain p p -[ ] , 3 with an equivalent one consisting of two square prisms D 1 and D 2 . In the former we assume